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This is a report of research conducted during the final period under NASA Grants 
NAGW-1937, "Tectonic History of the Terrestrial Planets". During this period the grant has 
supported the scientific contributions of a principal investigator, one Postdoctoral Associate 
(Suzanne E. Smrekar), and two Graduate Research Assistants (Patrick J. McGovern and 
Noriyuki Namiki) to the NASA Planetary Geology and Geophysics Program. 

The simplest way to convey the full scope and principal results of research during the 
past year is through the papers and abstracts of oral presentations supported by this grant 
during the final period. Following is a list of 5 papers, one book chapter, and 14 abstracts 
and extended abstracts that have been published or submitted for publication. This list serves 
as a table of contents for the bulk of this report. 


TABLE OF CONTENTS 


(PAPERS AND ABSTRACTS SUPPORTED BY NASA GRANT 
NAGW-1937 DURING APRIL 1992 - MARCH 1993) 


1. Evans, S.A., M. Simons, and S.C. Solomon, Flexural analysis of uplifted rift flanks on Venus, 

in Papers Presented to the International Colloquium on Venus , Pasadena, Calif., pp. 30-32, 
Lunar and Planetary Institute, Houston, Tex., 1992. 

2. McGovern, PJ., and S.C. Solomon, Patterns of deformation and volcanic flows associated 

with lithospheric loading by large volcanoes on Venus, in 1992 Spring Meeting, Eos Trans. 
Amer. Geophys. Un., 73, Suppl., 178, 1992. 

3. McGovern, P.J., and S.C. Solomon, Estimates of elastic plate thicknesses beneath large 

volcanoes on Venus, in Papers Presented to the International Colloquium on Venus, 
Pasadena, Calif., pp. 68-70, Lunar and Planetary Institute, Houston, Tex., 1992. 

4. McGovern, P.J., and S.C. Solomon, Aspects of modelling the tectonics of large volcanoes on 

the terrestrial planets, in Lunar and Planetary Science XXTV, pp. 959-960, Lunar and 
Planetary Institute, Houston, Tex., 1993. 

5. McGovern, P.J., and S.C. Solomon, State of stress, faulting, and eruption characteristics of 

large volcanoes on Mars, J. Geophys. Res., submitted, 1993. 

6. McKenzie, D., P.G. Ford, C. Johnson, D. Sandwell, B. Parsons, S. Saunders, and S.C. 

Solomon, Features on Venus generated by plate boundary processes, in 1992 Spring 
Meeting, Eos Trans. Amer. Geophys. Un., 73, Suppl., 304, 1992. 

7. Namiki, N., and S.C. Solomon, The gabbro-eclogite phase transition and the elevation of 

mountain belts on Venus, in Papers Presented to the International Colloquium on Venus, 
Pasadena, Calif., pp. 74-76, Lunar and Planetary Institute, Houston, Tex., 1992. 

8. N amiki , N. and S.C. Solomon, The gabbro - eclogite phase transition and the elevation of 

moun tain belts on Venus, J. Geophys. Res., submitted, 1992. 

9. Schubert, G., S.C. Solomon, D.L. Turcotte, M.J. Drake, and N.H. Sleep, Origin and thermal 

evolution of Mars, in Mars, edited by H.H. Kieffer, B.M. Jakosky, C.W. Snyder, and M.S. 
Matthews, Univ. Arizona Press, Tucson, pp. 147-183, 1992. 

10. Simons, M., B.H. Hager, and S.C. Solomon, Geoid, topography, and cgnvection-driven 

crustal deformation on Venus, in Papers Presented to the International Colloquium on 
Venus, Pasadena, Calif., pp. 1 10-112, Lunar and Planetary Institute, Houston, Tex., 1992. 

11. Simons, M., B.H. Hager, and S.C. Solomon, Geoid, topography, and convection-driven 

crustal deformation on Venus, in Lunar and Planetary Science XXIV, pp. 1307-1308, Lunar 
and Planetary Institute, Houston, Tex., 1993. 


12. Simons, M., and S.C. Solomon, Geoid-to-topography ratios on Venus: A global perspective, 

in 1993 Spring Meeting, Eos Trans. Amer. Geophys. Un., 74, Suppl., in press, 1993. 

13. Smrekar, S.E. and S.C. Solomon, Gravitational spreading of high terrain in Ishtar Terra, 

Venus,/. Geophys. Res., 97, 16,121-16,148, 1992. 

14. Smrekar, S.E., and S.C. Solomon, Constraints on crustal rheology and age of deformation 

from models of gravitational spreading in Ishtar Terra, Venus, in Papers Presented to the 
International Colloquium on Venus, Pasadena, Calif., pp. 1 14-116, Lunar and Planetary 
Institute, Houston, Tex., 1992. 

15. Solomon, S.C., Tectonic processes on Venus: Comparisons and contrasts with the Earth, 

Amer. Assoc. Petrol. Geol. 1992 Annual Convention Official Program, Calgary, Canada, p. 
124, 1992. 

16. Solomon, S.C., The tectonics of Venus: An overview, in Papers Presented to the 

International Colloquium on Venus, Pasadena, Calif., pp. 118-119, Lunar and Planetary 
Institute, Houston, Tex., 1992. 

17. Solomon, S.C., The tectonic and volcanic evolution of Venus: Catastrophic or gradual?, in 

1992 Fall Meeting, Eos Trans. Amer. Geophys. Un., 73, Suppl., 328-329, 1992. 

18. Solomon, S.C., Venus: Keeping that youthful look. Nature, 361, 1 14-1 15, 1993. 

19. Solomon, S.C., A tectonic resurfacing model for Venus, in Lunar and Planetary Science 

XXTV, pp. 1331-1332, Lunar and Planetary Institute, Houston, Tex., 1993. 

20. Solomon, S.C., The resurfacing controversy for Venus: An overview and a mechanistic 

perspective, in 1993 Spring Meeting, Eos Trans. Amer. Geophys. Un., 74, in press, 1993. 

21. Squyres, S.W., D.G. Jankowski, M. Simons, S.C. Solomon, B.H. Hager, and G.E. McGill, 

Plains tectonism on Venus: The deformation belts of Lavinia Planitia, J. Geophys. Res., 
97, 13,579-13,599, 1992. 


FLEXURAL ANALYSIS OF UPLIFTED RIFT FLANKS ON VENUS; Susan A. 
Evans, Mark Simons, and Sean C. Solomon, Department °f Ea r th, _Atm°sP h ' nc ' and Planetary 
Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139. 


Introduction. Knowledge of the thermal structure of a planet is vital to a thorough 
understanding of its general scheme of tectonics. Since no ^ect meas^em^tsofheat 
thermal gradient are available for Venus, most estimates have been derived fr ® m “eoret cal 
considerations or by analogy with the Earth [1]. The flexural response of the ^osphereto apphed 
loads is sensitive to regional thermal structure. Under the gumption i that t the yield strcng 
function of depth can be specified, the temperature gradient can be inferred from the efi Fecttye 
elastic plate thickness [2]. Previous estimates of the effective elastic plate ^ k " eS ^ n J C ^ 
range from 11-18 km for the foredeep north of Uorsar Rupes [3] to 30-60 km lor the wmular 
troughs around several coronae [4,5]. Thermal gradients inferred for these regions are 14-23 K 
km-1 and 4-9 K km-l , respectively [3,4]. In this study, we apply the same techniques to 
investigate the uplifted flanks of an extensional rift. Hypotheses for the on ^ n ^ ft ® d .^ a f 
flanks on Earth include lateral transport of heat from the center of the nft, vemcal a^s^ of heat 
by small-scale convection, differential thinning of the lithosphere, dynarrucal uplift, a^d ^static 
response to mechanical unloading of the lithosphere [6]. The last hypothesis is considered the 
dominant contributor to terrestrial rift flanks lacking evidence for volcanic 
rift structures which are no longer active [6]. In this study, we model the uplifted flanks of a 

Venusian rift as the flexural response to a vertical end load. in an area to 

Tectonic Environment. We examine a linear nft system centered at 33 S, 92 E, in an area to 
the east of Aino Planitia. The feature appears as a linear ridge in Pioneer Venus 
has been named Juno Dorsum. However, the increased resolution of Magellan images and 
topography has established that this feature is actually a linear rift with bighs^ 

The rift h 100 km wide and 450 km long; it has a central depression 1-2 km deep ^d flanks 
elevated by as much as 1 km (Figures 1 and 2). The rift connects a 1.5-km-high volcano on its 
western edge to two coronae, Tai Shan and Gefjun [7], to the east Despite the presence of the 
volcano Sd coronae, both the center of the rift and the adjacent flanks appear to be free ofvocan^c 
flows. Juno Dorsum is at the end of a nearly continuous chain of coronae, nfts, and linear 
fractures which extends eastward and appears to terminate at the northwestern edge of Artemis 
Corona. The rift appears in Magellan images as numerous ENE-trendmg lineaments which we 
interpret to be normal faults. To the north and south are smooth, radar-dark plains which stand 
very close to mean planetary elevation. 
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Figure 1. Topographic contour map of the 
Juno Dorsum region. Thin solid lines denote 
positive elevation contours, thick solid lines 
zero elevation, and dashed lines negative 
elevations. The contour interval is 0.3 km, 
and the datum is mean planetary radius, 
6051.9 km. The rift is located in the center 
right of the figure, with an unnamed volcano 
to its west. The north-south-trending lines 
indicate the tracks of orbits 965-968 and 971- 
972 (numbers increasing eastward). The 
orbit tracks make an angle of approximately 
75* with the strike of the rift 
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Figure 2. Example orbit profiles over Juno Dorsum; see Figure 1 for location. The rift flanks are 
elevated 800 to 1000 m above the surrounding plains, 

topographic features is minimal. Orbit 967 has been offset vertically by 2500 m for clanty. 

We use Magellan altimetric profiles from orbits 965-968 and 970-972 figures 1 2). 

While topographic profiles perpendicular to the strike of the rift are preferable, the orbitffacks 
cross the riffatan angle of approximately 75*, and we can easily make the geomemcal comecuon 
for this small difference. Profiles from single orbits avoid errors introduced by inter-orbit o s . 
As is evident in Figure 2, the topography of the surrounding plains region is rather smooth, and 
Ae uplift of die flanks is quite distinct Li easy to discern from other regional topography. 

Model We model the uplift of the rift flanks as the flexural response of a thin elastic plate 
overlying a fluid substratum [8]. We assume that the flexure is of the form 

w = exp(j)[Acos(^] + Bsin(i)] (1) 

where w is the vertical deflection of the plate from the regional datum, x is the distance along the 
profile, and a is the flexural parameter given by 


a = 


Ell 


|_3pmg(l-v 1 2 )] 


( 2 ) 


where E is Young's Modulus (60 GPa), T c is the elastic plate thickness, p m is the > mantle d lensity 
(3300 kg nr 3) g is the gravitational acceleration (8.87 ms), and v is Poisson ( 
obiective is to find the tiiree parameters A, B, and a such that the root mean square (rms) misfit 
between the observed and calculated topographic profiles is a minimum. Since ^depends non- 
linrarlv on a we first fix the value of a, then calculate the best-fitnng values of A and B. The 
flexural parameter is then varied over a range corresponding to an elastic plate thickness of 1 to 50 
km. and die procedure is repeated. From the parameters providing the best overall fit we calculate 
the* plate ouvature at the ftet zero-crossing. This curvature and the apparent elastic plate thickness 

1 “0 «• We also remove the mean elevation 

and thelLt-fit linear trend from our profiles before fitting to prevent bias from long- wavelength 
regional topography. 
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Results. Of the fourteen profiles obtained in orbits 965-968 and 971-972, we andyze ten 
which have flanks rising 600-1200 m above the regional mean elevation. Four of *«e profiles 
are well fit by elastic plate thicknesses of 8-16 km. For each of these profiles, there is a narrow 
range of plate thicknesses for which the rms error is small, so the results arewell constrained 
(Figure 3). Three of these profiles produce best-fitting plate thicknesses of 8-9 km, suggesting that 
the true effective elastic plate thickness may be near the low end of this range. The remaining 
profiles generally contain more small-scale topographic variation of ± 100 m in amplitude and are 
equally well fit by a range of plate thicknesses. The fits for all profiles, however, are consistent 

with elastic plate thicknesses of 8-20 km. 
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Figure 3. Left: topographic profiles (solid) and best-fit flexural profiles (dashed) for three half 
profiles; top to bottom are orbit 967, southern flank; orbit 968, southern flank; orbit 968, northern 
flank. Right: the rms error as a function of elastic plate thickness for the same profiles. Note tha 
the top two profiles are best fit only over a small range of elastic plate thicknesses, whereas the 
bottom profile can be fit by a broader range of plate thicknesses. All profiles are consistent with 
elastic plate thicknesses in the range 8-20 km. For clarity, the upper two profiles on the left have 
been vertically offset by 1000 m and 2000 m, respectively, and the corresponding error curves by 

30 and 60 m. 

These effective plate thicknesses and the calculated curvatures imply thermal gradients of 15- 
30 K km* 1 if the base of the mechanical lithosphere is limited by the ductile strength of olivme 13J. 
Because the values of elastic lithosphere thickness derived here are comparable to or less than the 
estimated crustal thickness beneath plains regions on Venus [9,10], the base of the mechanical 
lithosphere may lie within the crust and be limited by the flow law for lower crustal material, 
such a situation, the thermal gradients would be less than the figures given above. 

The effective elastic plate thicknesses we denve for the Juno rift strocture (8-20 km) lare 
comparable to those obtained for the North Polar Plains at Uorsar Rupes (11-18 km) and both 
regions involve flexure of the lithosphere beneath plains at elevations near the planetary mean. A 
general relationship between regional elevation and elastic plate thickness might be expected if 
regional elevation were controlled by lithospheric thermal structure [11], as in the Earth s oceanic 
regions. Further analysis of flexure is warranted to determine if such a relationship holds tor 

Venus. 
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Patterns of Deformation and Volcanic Flows Associated 
with Lithospheric Loading by Large Volcanoes on Venus 

Patrick J. McGovern and Sean f- Sol^on (Dept of Earth, 
Atmospheric, and Planetary Sciences, Cambridge, MA 02139) 

Magellan radar imaging and topography data arc now avmlable for a 
number of volcanoes on Venus greater than 100 km m radius. These 
data can be examined to reveal evidence of the flexural response of the 
lithosphere to the volcanic load. On Venus, erosion and sediment 
deposition are negligible, so tectonic evidence of deformation around 
large volcanoes should be evident except where buried by very young 
flows. Radar images reveal that most tectonic features and flow units 
on the flanks of these volcanoes have predominantly radial orientations. 
However, both Tepev Mons in Bell Regio and Sapas Mons in Atla 
Regio exhibit circumferential graben on their flanks. In addition, 
images reveal several flow units with an annular character around the 
north and west flanks of Tepev Mons. This pattern most likely results 
from ponding of flows in an annular flexural moat Maat Mons in Atla 
Regio and Sif Mons in Eistla Regio are examples of volcanoes that lack 
circumferential graben and annular flows; discernible flow units and 
fractures on these constructs appear to be predominantly radial. 
Altimetry data can also provide evidence of flexural response. Tepev 
Mons is partially encircled by depressions that may be sections of a 
flexural moat that has not been completely filled. The lwatioM of these 
depressions generally coincide with the annular flows described above 
There is weaker evidence for such depressions around Maat Mons as 
well. The lack of circumferential tectonic features around most 
volcanoes on Venus might be explained by gradual moat filling and 
coverage by radial flows. The depressions around Tepev (and possibly 
Maat) may indicate that this process is currently continuing. We use 
analytic models of plate flexure in an axisymmetnc geometry to 
constrain the elastic plate thickness supporting Tepev Mons. If we 
consider the outer radius of the ponded flows to be the edge of a moat, 
we find that models with elastic plate thickness of 10-20 km fit best. 
Finite element models of a volcanic load detached from the underlying 
lithosphere predict overthrusting and radial normal faulting at the 
volcano's edge. Such a mechanism for the formation of radial nft zones 
on Venus volcanoes would make such features analogous to structures 
on the flanks of volcanoes on Earth. 


ESTIMATES OF ELASTIC PLATE THICKNESSES BENEATH LARGE 
VOLCANOES ON VENUS; Patrick J. McGovern and Sean C. Solomon, Dept, ot fcartn, 
Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 
02139. 

Introduction. Magellan radar imaging and topography data are now available for a 
number of volcanoes on Venus greater than 100 km in radius. These data 
reveal evidence of the flexural response of the lithosphere to the volcanic load. On Ear^flexure 
beneath large hotspot volcanoes results in an annular topographic moat which is paroally to 
completely filled in by sedimentation and mass wasting from the volcano s flanks (see [1, 2]). Un 
Venus, erosion and sediment deposition are considered to be negligible at t he r esolution of 
Magellan images [3]. Thus, it may be possible to observe evidence of flexure by the ponding of 
recent volcanic flows in the moat. We also might expect to find topographic signals fromu^iUed 
moats surrounding large volcanoes on Venus, although these signals may be partially o^scuredby 
regional topography. Also, in the absence of sedimentation, tectonic evidence of deformatio 
around large volcanoes should be evident except where buried by very young flows. 

We have found two examples to date of volcanoes with strong evidence for moat formation 
and infilling by flows. Radar images of Tepev Mons, a volcano about 125 km m radius .near the 
southwestern comer of Bell Regio, reveal a bright flow unit draped around the northern and 
SSb ofthe volcano. An Unnamed volcano at 10°N, 275°E (southwest of Beta Regio) 
exhibits both circumferential flanking flows and lateral spreading of originally radially trending 
flows The distal edges of these flows terminate about 240-260 km from the summit The edges of 
these flows form an arc of greater than 120° to the north and west of the volcano. 

Method. We use analytic solutions in axisymmetric geometry [4] for deflections and stresses 
resulting from loading of a plate overlying an inviscid fluid. Solutions for a set of disk loads are 
superimposed to obtain a solution for a conical volcano. The deflection of the hthosphere produces 
an annular depression or moat, the extent of which can be estimated i by i measuring I ^^stance 
from the volcano's edge to the first zero crossing or to the peak of the flexural inch. Magellan 
altimetry data records (ARCDRs) from data cycle 1 are process^ using i^e GMT mapping and 
graphics software [5] to produce topographic contour maps of the ta ^ 

topographic profiles that cut across the annular and ponded flows seenon the § , ‘ w 

comparing the locations of these flows to the predicted moat locations from a range of models, we 
wtiinate the elastic plate thickness which best fits the observations, together with the uncertainty in 

that estimate. 

Results. Figure 1 shows two cross sections through Tepev Mons. The areas covered by 
annular flows are marked with arrows. Figure 2 shows deflections calculated analyticaUy for a 
conical load of height 10 km and radius 150 km for elastic plate thicknesses T e of 10km « and 20 
km. Arrows denote the predicted approximate extent of a moat due to loading in each model. Note 
that for the analytic solutions, increasing the elastic plate thickness increases both the maximum 
depth and the radial extent of the predicted moat. The model in Figure 2a matches the no^western 
cross section in Figure la. The area of flows in Figure lb is somewhat larger in radial extent and 
is better matched tiythe predicted moat in Figure 2b. Figure 3 shows a cross section through the 
volcano at 10°N 275°E. Figure 4 shows calculated deflections for a conical load of height 5 km 
and radius 250 km, for elastic plate thicknesses of 10 and 20 km. The extent of annular flows is 
better matched by the smaller thickness. 

Discussion These two volcanoes are apparently atypical of large volcanoes on Venus in that 
both topography and flow morphology suggest the existence of a flexural moat. Onotherlarge 
shield volcanoes such as Sif Mons and Sapas Mons, topographic evidence of a moat is lada g, 
and no flows that have convincingly ponded in annular moats can be identified. Large volcanoes 
Aat fonn on or near large rift zones (such as Maat Mons, Ozza Mons, and others) also lack 
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topographic evidence of a flexural moat It should be noted that as altimetry data and stereo 
imaging from later cycles become available, coverage will improve, gaps will be filled in, and it 
may be possible to identify topographic signatures of moats which have eluded our search to date. 

In contrast to Tepev Mons and the construct at 10°N, 275°E, most large volcanoes on Venus 
are generally characterized by fractures and flow units that have dominandy radial onenta ons. 
Given this observation and the assumption that lithospheric flexure has occurred, then the moats 
must be filled or covered. Mass wasting, through large landslides or slumps such as those 
discovered off the Hawaiian Islands [2], is one possibility. Moat-filling flows may be covered by 
later radial flows, a process that may have been repeated many tunes during volcano growth. 

Conclusions. Under the assumption that circumferentially ponded flows surrounding 
volcanoes have at least partially filled a flexural moat, we can compare the posmon^and extent of 
the moat to predictions from analytic flexure models. For Tepev Mons, models tiiat fit best ^ha e 
elastic plate thicknesses in the range of 10 to 20 km. For the volcano at 10 N, 275 E, die best fit is 
obtained for a thickness of about 10 km. These values of elastic plate tinckness can be converted 
to estimates of mechanical plate thickness and thermal gradient [6]. Using parameters appropriate 
for Venus [7] we obtain thermal gradients values of about 12-25 K/km for Tepev Mons and about 
25 K/km for 10°N, 275°E. These gradients are in the range expected if Venus loses most of its 
internal heat by conduction through a globally continuous, if laterally heterogeneous, lithospheric 

shell [7]- 
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270, 1991; [4] J. F. Brotchie, Mod. Geol., 5, 15, 1971; [5] P. Wessel and W. R F. Smi , o 
Trans . AGU, 72, 441, 1991; [6] M. K. McNutt, J. Geophys. Res., 89, 11,180, 1984, [7] S. C. 
Solomon and J. W* Head, Geophys . Res. Lett., 17, 1393, 1990* 



Figure 1. Topographic cross sections through Tepev Mons. (a) Northwestern profile; (b) 
western profile. 
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Figure 3. Topographic cross section 


through the volcano at 10 “N, 275°E, extending westward. 
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ASPECTS OF MODELLING THE TECTONICS OF LARGE VOLCANOES ON 
THE TERRESTRIAL PLANETS; Patrick J. McGovern, Dept, of Atrno s phenc, and 
Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, and bean U. 
Solomon, Dept, of Terrestiial Magnetism, Carnegie Institution of Washington, Washington, D. v_. 

20015. 

Introduction. Analytic solutions for the response of planetary lithospheres to volcanic loads 
have been used to model faulting and infer elastic plate thicknesses [e.g., 1-4]. Predictions of the 
distribution of faulting around volcanic loads, based on the application of Anderson s cntena for 
faulting [5] to the results of the models, do not agree well with observations [6,7], buch models 
do not give the stress state in the load itself, but only suggest a state of horizontal compressive 
stress there. Further, these models have considered only the effect of an instantaneously emplaced 
load. They do not address the time evolution of stresses, nor do they consider the effect of a load 
which grows. A finite element approach allows us to assign elements to the load itself, and thus 
permits calculation of the stress state and stress history within the edifice. The effects of episodic 
load growth can also be treated. When these effects are included, models give much better 

Method. We use the finite element code TECTON [8,9] to construct axisymmetnc niodels of 
volcanoes resting on an elastic lithospheric plate overlying a viscoelastic asthenosphere. We have 
implemented time-dependent material properties in order to simulate incremental volcano growth. 
The viscoelastic layer was taken to extend to a sufficient depth so that a ngid lower boundary has 
no significant influence on the results. The code first calculates elastic deformations and stresses 
and then determines the time-dependent viscous deformations and stresses. Time in the model 
scales as the Maxwell time x m in the asthenosphere. ... , 

Results We consider a volcano 25 km in height and 200 km in radius on an elastic 
lithosphere 40 km thick (parameters approximately appropriate to Ascraeus Mons). The volcano 
consists of three load increments applied at intervals of 1000 Xm- Contours of maximum deviatonc 
stress in the fully-grown edifice at die conclusion of flexure (t = 3000 x m ) are shown in Figure 1. 
Note that the maximum stress occurs in the lower mid-section of the edifice. We adopt the 

convention that tension is positive. . . , t 

Discussion. Simple models of plate flexure due to an instantaneous load predict a state of 
horizontal compressive stress in the plate beneath the load [1,4] with an increasing magnitude of 
compression from the middle of the plate upward, and a maximum at the surface. Models witha 
succession of loads emplaced from the bottom up display a different distribution of stress. 1 he 
maximum horizontal stress occurs in the lowest (first) load increment and decreases upward with 
each later load increment. Thus, thrust faults (predicted by simpler mode s) associated with such 
. stresses would be expected to occur only early in the evolution (or not at all) and would be covered 
by later units which remain unfaulted. Such a distribution of stresses could also affect the 
locations and dimensions of magma chambers. A zone of horizontal compression within the 
edifice would inhibit ascending magma from reaching the summit and thus might cause magma to 
pool beneath it. Since magma will propagate perpendicular to the direction of least compressive 
stress, radial propagation of magma in sheet dikes or sills might also occur. The effects of adding 
fault slip along a basal detachment between edifice and lithosphere (as proposed in models of the 
flank tectonics of Hawaii [10]) was shown to achieve an even greater reduction in horizontal 
compressive stress in the edifice [11]. Such a structure may control the growth and evolution of 
tectonic features surrounding volcanoes on both Mars and Earth [11]. t . 

Stresses calculated from planetary loading models have been used, in conjunction with 
Anderson’s criteria for faulting [5], to predict types of faulting that should be observed [1,2,4] 
Typically a zone is found surrounding the load with both the least compressive stress Oi and the 
most compressive stress o 3 horizontal (assuming that near the surface, o« is small and 
compressive). Application of Anderson's criteria yields a prediction of stnke-slip faulting. Such a 
region surrounding a large load has not heretofore been observed [6,7]. Given the complete stress 
field in an element in this region from the TECTON solutions, we can plot Mohr diagrams, as in 
Figure 2, and compare them with a Mohr-Coulomb failure envelope. The Mohr circles grow in 
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size as flexure proceeds, until the failure envelope is exceeded at about 70 tm ^ tcr 
The largest circle shows stresses at the conclusion of flexure. This stress suite <^^sponds to that 
given by analytic plate models (in which flexure occurs instantaneously). Sunply int<^retog dus 
s«ss stare is P seen to be misleading for two reasons. First, 9 , i<« *«) * l^e^ die 

tensional strength of any rock. Use of a shear failure criterion (such ^ Mo^-Ojdomb) for such a 
stress state is improper [7]. Use of a tensile failure criterion will. predict the formation of tensmn 
joints which may be precursors of graben [7]. Second, consideration of the stress history mdicate 
that at the time of first failure, the shear failure criterion is satisfied, with principal s «ss dumions 
that predict circumferential normal faulting. Formation of these faults wou . .j ^ 

provide planes of low resistance along which further faulting with a similar 2 

expected to develop. Thus, the "prediction" of a zone of ^ke-slip faulting^ sirrrounding a 
lithospheric load, based on simple instantaneously-loaded plate flexure models, 1S m eiror. 

Tonc£bns. Simple interpretations of the results ^m mstantaneous load modds lead to 
predictions of lithospheric stress fields that do not match observations Constder^onof^e 
history of stresses in the lithosphere and within the load itself helps to resolve these discrete . 
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Figure 1. Contours of maximum deviatoric stress (MPa) in the volcanic edifice described above. 
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Figure 2 Mohr diagram for element at the top of the lithosphere (stresses calculated at element 
S? r = 265 km T= -2.5 km). Diagonal tine is Mohr-Coulomb failure envelope for basalt [12] 
Srt. amde of internal friction <t> = 49° and cohesion c = 38 MPa. In order of increasing radius, the 
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first exceeded (r = 70 t m ), Gi = c„ and a 3 =o K . At final time shown, c x - c„ and o 3 - Gee- 
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ABSTRACT 

The formation of a large volcano loads the underlying lithospheric plate and can lead to 
lithospheric flexure and faulting. In turn, lithospheric stresses affect the stress field beneath and 
within the volcanic edifice and can influence magma transport Modeling the interaction of these 
processes is crucial to an understanding of the history of eruption characteristics and tectonic 
deformation of large volcanoes. We develop models of time-dependent stress and deformation for 
the Th arsis volcanoes on Mars. A finite element code is used that simulates viscoelastic flow in the 
mant le and elastic plate flexural behavior. We calculate stresses and displacements due to a 
volcano-shaped load emplaced on an elastic plate. Models variously incorporate growth of the 
volcanic load with time and a detachment between volcano and lithosphere. The models illustrate 
the manner in which time-dependent stresses induced by lithospheric plate flexure beneath the 
volcanic load may affect eruption histories, and the derived stress fields can be related to tectonic 
features on and surrounding Martian volcanoes. After an initial load increment is emplaced, 
flexurally-induced stresses grow with time and the principal stress directions in the volcano rotate 
as flexure proceeds. Magma is expected to propagate perpendicular to the least compressive stress 
axis. As a result of flexure, this axis rotates from a horizontal orientation to a nearly vertical one, 
thus magma propagation paths will tend to rotate from vertical to horizontal orientations. We 
suggest that at the later stages of flexure, this effect would tend to favor eruption sites on to the 
flanks of the volcano rather than the summit. Such a scenario is consistent with the 
photogeologically determined evolution of the Tharsis Montes. As a result of flexure there are 
three regions where stresses become sufficiently large to cause failure by faulting (according to the 
Mohr-Coulomb criterion): at the surface of the plate just outward of the volcano, and at the base of 
the elastic lithosphere beneath the center of the volcano, and on the upper flanks of the first 
volcanic load increment. Normal faulting is the dominant mode of failure predicted for the first 
region, consistent with circumferential graben observed around the Tharsis Montes and with the 
scarp at the base of Olympus Mons, interpreted as a large-offset, listric normal fault. Normal 
faulting, mostly radially oriented, is predicted for the second region. Under the premise that failure 
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in the sub-volcano zone may be strongly influenced by pre-existing weak zones or regional 
stresses, this feature may have a surface expression in the rifting and the bilateral symmetry of the 
Tharsis Montes. Failure in the third region is predicted to consist of thrust faulting, oriented 
mostly circumferentially on the upper flanks with an annulus of radially oriented thrust faulting 
starting about midway downslope. Concentric terraces, interpreted by some workers as thrust 
faults, on the upper flanks of Olympus Mons may correspond to the predicted circumferential 
thrust features. For volcanoes detached from the plate, predicted failure in the edifice takes the 
form of radial normal faulting at the base. Addition of an extensional stress due to regional 
topography yields a pattern of predicted faulting which closely matches that- observed on the 
Tharsis Montes. This stress state is also consistent with an interpretation of the aureole deposits of 
Olympus Mons as the result of gravity sliding along a basal detachment. 

Our models also suggest explanations for the lack of strike-slip features, predicted by simple 
flexural models, around the Tharsis volcanoes. For a given load increment, the first mode of near- 
surface failure for most of the area immediately outward of the load is circumferential normal 
faulting and graben formation. As the volcano grows and the flexural response to the increasing 
load proceeds, the predicted failure mode in a portion of this annular region surrounding the 
volcano changes to strike-slip faulting. Because normal faulting has been predicted to have taken 
place earlier, however, it is likely that release of later stresses will occur by reactivation and 
growth of these normal faults and graben rather than by the formation of strike-slip faults. 
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INTRODUCTION 

The major Tharsis volcanoes of Mars (Olympus Mons and the three Tharsis Montes. Arsia 
Mons, Pavonis Mons, and Ascraeus Mons) are among the largest known volcanic structures in the 
solar system. An understanding of the formation and evolution of these structures can provide 
important constraints on the processes that built and maintained the Tharsis Rise. In addition, the 
Tharsis volcanoes may be analogues to large hot-spot volcanoes on Earth, such as those of 
Hawaii. For instance, Kilauea and Olympus Mons have very similar ratios of volcano height to 
diameter and of basal scarp height to volcano height [Borgia et al., 1990]. Thus studies of the 
evolution of large Martian volcanoes may yield insight into terrestrial volcanic processes as well. 

In this paper we utilize finite element models to evaluate the evolution in internal stress and 
deformation within and surrounding the Tharsis volcanoes, and we discuss how the time- 
dependent stress field may be relatable to the eruption characteristics of the volcanoes and to the 
formation of associated tectonic features. 

To date the investigation of the evolution of stresses in large volcanoes has taken two paths, 
models of edifice stresses alone (usually finite element models with rigid bottom boundary 
conditions), and investigations of flexural stresses in the lithospheric plate supporting the volcano. 
As examples of the first category of study, Chevallier and Verwoerd [1988] used an axisymmetric 
planar finite element code to investigate the effect of magma chamber and external pressures on 
stress and eruption histories of hot spot volcanoes, Dieterich [1988] modeled stress in volcano rift 
zones by means of a two-dimensional triangular grid of elements, and Ryan [1988] employed a 
horizontal planar finite element model of the flank of Kilauea to determine displacements and 
stresses due to dike emplacement. In an example of the second class of study, Thurber and Gripp 
[1988] applied a flexural model to constrain the tectonics of volcano flank motions, ten Brink and 
Brocher [1987] proposed a link between flexural stresses in the lithosphere and eruption history; in 
their scenario the orientation of flexural stresses beneath a given point on the volcanic chain 
changes with time as the volcanic load is emplaced and then eroded. Given that magma seeks to 
propagate along paths perpendicular to the direction of least compressive stress, magma ascent can 
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be blocked when both principal horizontal stress deviators are compiessional. Thus eruption 
history depends on location along the volcanic chain or on time since volcano formation. 

Recent studies have also been carried out on edifice stresses of large Martian volcanoes. 
Thomas et al. [1990] investigated the tectonics of the flanks of Martian volcanoes by means of an 
incompressible finite-element model. Volcano self-loading and magma chamber effects were 
included, but lithospheric flexure was not. Thomas and coworkers found that stresses on the 
flanks of a large volcano are sufficient to cause circumferentially oriented thrust faulting, and they 
suggested that such thrusting produced the concentric terraces observable high on the flanks of 
Olympus Mons. Zuber and Mouginis-Mark [ 1992 ] treated the evolution of the surface of 
Olympus Mons, using a finite element model to calculate stresses in the summit caldera region 
caused by different magma chamber locations and geometries for comparison with observed 
patterns of faulting. 

In this paper we study the stress field within a volcano and the lithosphere upon which it rests 
as a unified system. With such a model formulation we can account explicitly for the interaction 
between edifice stresses (e.g., due to volcano self-loading) and flexural stresses (a result of the 
load induced by the volcano on the lithosphere). We include both growth of the volcano and 
viscoelastic deformation in the asthenosphere, so the problem is intrinsically time-dependent. The 
calculated displacements yield the subsidence history of the volcano. The orientations of principal 
stresses and their change with time provide important constraints on possible magma emplacement 
paths and eruption histories. With the computed stress fields, a failure criterion can be used to 
predict locations and modes of faulting within and near the volcano. After a short discussion of 
important geological and geophysical characteristics of the Tharsis volcanoes, we briefly describe 
the finite element procedure and the modelling assumptions. The results of the numerical 
computations are next presented, and their potential implications for the evolution of eruption 
characteristics and the formation of tectonic features are then compared with known constraints on 
the evolution of the Tharsis volcanoes. 
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CHARACTERISTICS OF MAJOR THARSIS VOLCANOES 

The four major Tharsis volcanoes are the largest of the martian shield volcanoes [e.g., Greeley 
and Spudis, 1981]. Each construct is composed of many overlapping flows and flow units 
erupted over a period of activity as much as 2-3 Gy. in duration [Tanaka, 1986]. The approximate 

heights and basal diameters are given in Table 1. 

The largest of the major Tharsis volcanoes is Olympus Mons (Figure 1). The principal 
tectonic features associated with Olympus Mons are (1) the summit caldera complex, consisting of 
a series of circular depressions with complex patterns of faulting, (2) concentric terraces seen on 
the upper slopes of the volcano, (3) the basal scarp, a nearly vertical cliff surrounding the volcano 
at a radius of about 300 km, and (4) the aureole deposits, which occupy a vast region dominantly 
to the northwest of the main shield (downslope from the Tharsis rise [U. S. Geological Survey 

1989] ). The summit caldera is likely the result of a collapse following withdrawal of magma from 
a high-level magma chamber [Mouginis-Mark, 1981; Zuber and Mouginis-Mark, 1992]. As noted 
earlier Thomas et al. [1990] have interpreted the concentric terraces to be thrust faults. The basal 
scarp has been variously interpreted as a thrust fault [Morris, 1981], a listric normal fault [Francis 
and Wadgt, 1983], and a fault-propagation fold over a subsurface thrust fault [Borgia et al., 

1990] . The aureole deposits are generally held to be disrupted landslide material derived from the 
slopes of the volcano [Harris, 1977; Lopes et al, 1980; Francis and Wadge, 1983]. Tanaka 
[1985] proposed that gravity sliding and spreading of the aureole deposits is enabled by a weak 
basal detachment between aureole and substrate. A layer containing about 10% interstitial or 
interbedded ice would be sufficiently weak to provide such a detachment. 

Tectonic features on the flanks of Hawaiian volcanoes on Earth may provide useful analogues 
to the Olympus Mons scarp and aureole. Acoustic backscatter images of the Hawaiian Ridge 
reveal that extensive mass wasting deposits (slumps and debris avalanches) flank all of the major 
Hawaiian shields [Moore et al., 1989]. A study of focal mechanisms of shallow earthquakes on 
the southern flanks of Kilauea [Thurber and Gripp, 1988] supports the hypothesis that a 
detachment surface separates the volcanic edifice and the older oceanic crust. By this hypothesis, 
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the volcano can overthrust the underlying crustal layer, sliding on a decollement of weak oceanic 
sediments. The slides and slumps, driven by intrusions at the rift zones, detach from the volcano 

along large-offset listric normal faults [Lipman et al., 1985]. 

Tectonic features observed on and around the Tharsis Montes volcanoes differ somewhat from 
those of Olympus Mons. The principal tectonic features associated with these volcanoes are (1) 
circumferential graben on the lower volcano slopes and the surrounding plains, (2) linear 
embayments or rifts approximately bisecting each volcano along an axis trending N40°E and 
serving as sources for flows embaying the northeast and southwest flanks, and (3) lobe shaped 
deposits to the west or northwest of the edifice [Scott and Tanaka, 1981; Scott et al, 1981a, b, c]. 
The Tharsis Montes volcanoes exhibit an approximate bilateral symmetry about a NE-SW-trending 
axis coinciding approximately with the line connecting their centers. Graben tend to occur on the 
northwest and southeast flanks of these volcanoes (Figures 2-5), i.e. perpendicular to radial 
directions broadly orthogonal to the line of the rifts. From photogeological study of Viking Orbiter 
images, Crumpler andAubele [1978] proposed the following evolutionary sequence for the 
Tharsis Montes: (1) construction of the main shield, (2) outbreak of parasitic eruption centers on 
the volcano flanks along the NE-SW-trending axis, (3) subsidence of the summit and formation of 
concentric fractures and graben, and (4) formation of a bisecting rift along the NE-SW-trending 
axis, with rift eruptions leading to flooding of the summit depression and inundation of the rifted 
flanks. This sequence is most advanced on Arsia Mons; Pavonis Mons has reached stage (3), and 
Ascraeus Mons stage (2). The morphology and tectonics of the flanks of the Tharsis Montes are 
described in greater detail by Zimbleman and Edgett [1991]. They propose that the lobate deposits 
formed by gravity sliding and were further modified by both effusive and pyroclastic volcanic 

activity. 

Prominent rift zones that radiate from the summit calderas are dominant structural features 
characteristic of Hawaiian shield volcanoes [Peterson and Moore, 1987]. These zones are sites of 
emplacement of magma drained from the summit magma chamber and indicate an environment of 
horizontal extension. The rift zones and landslides or slump deposits appear to be causally related. 
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as described above. The Hawaiian rift zones may be analogues to the radial rifts of the Tharsis 
Montes. 

The thickness of the elastic lithosphere beneath the Tharsis Montes volcanoes has been 
inferred from the radial distances of their circumferential graben by Comer et al. [1985] (see 
Figures 3-5). Preferred values of elastic lithosphere thickness are around 20 km. Concentric 
graben are not found around Olympus Mons, however, which led Comer and coworkers to 
conclude that the thickness of the elastic lithosphere beneath Olympus Mons must be much greater 
(> 150 km) than beneath the Tharsis Montes. Of course, the actual lithosphere does not behave 
perfectly elastically; rather its strength is limited by frictional failure at shallow depth and by ductile 
flow at greater depth. From strength envelope considerations for crustal and mantle material 
[McNutt, 1984], values for elastic plate thickness T e were converted to estimates of mechanical 
plate thickness T m and lithospheric thermal gradient by Solomon and Head [1990]. For Olympus 
Mons, flexurally induced curvature is small and T m is approximately equal to T e . The lithospheric 
thermal gradient is then less than 5 K km-*. For the Tharsis Montes, flexural curvatures are larger, 
T m exceeds T e , and mean thermal gradients (corresponding to the best fitting values of T m ) are in 
the range 10-14 K km' 1 . Upper and lower bounds on T c for the Tharsis Montes from Comer et al. 
[1985] allow a range of thermal gradients from 7-27 K km* 1 . 

METHOD 

We use the finite element program TECTON [Melosh and Rafesky, 1980, 1983] to model 
stresses and displacements in a large volcano and in the crust and mantle beneath and around the 
volcano. TECTON's capability for modeling a viscoelastic rheology in the mantle allows us to 
model the time-dependent flexural response of an elastic lithosphere to each major increment of 
load. The program first calculates the elastic (i.e., "instantaneous") response to a load. The 
stresses and displacements arising from load-induced viscoelastic flow in the mantle are then 
calculated for a specified number of time steps. The Maxwell time (x M ) of the mantle, defined as 
the ratio of viscosity ri to shear modulus |i, is used as a convenient reference time scale (in 
subsequent discussion the term "Maxwell time" will refer to the mantle Maxwell time). Complete 
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descriptions of stresses and displacements were obtained at 5, 10, 15, 20, 30, 40, 70, 100, 200, 
400, 700 and 1000 Maxwell times after emplacement of each load increment. 

We assume that the problem is axisymmetric, with cylindrical coordinates r, 0, and z. Out-of- 
plane shear stresses and o 0z are then zero. We solve for r and z displacements and stress 
components a m a 0e , a^, and a n . In axial symmetry, two principal stresses are confined to the rz 
plane, and a 00 , the stress normal to this plane, is also a principal stress. An example of the finite 
element grid used for this study is shown in Figure 6. The displacement boundary conditions are 
that nodes on the side walls (r = 0 and r = r max ) are fixed in r but free to move in z, and that nodes 
on the bottom boundary (z = -820 km in the example shown) are fixed in z but free to move in r. 
The lower comers of the box are fixed in both directions. The volcano in the example has a radius 
of 200 km and is 25 km in height, the approximate present dimensions of Ascraeus Mons. In 
Figure 6 the volcano is the triangular region in the upper left-hand comer. The volcano rests on 
top of an elastic lithospheric plate of thickness T e . All elements in the volcanic edifice and the plate 
have a high viscosity appropriate to the lithosphere. These elements behave essentially elastically 
over the time scales considered here. All elements below the elastic lithosphere have a lower 
viscosity value, appropriate for an asthenosphere. These elements experience viscous relaxation 
over the time modeled. Material property values adopted in our calculations are listed in Table 2. 
Parameters such as density and Young's modulus differ for the crust and mantle. Elements above 
depth tf. are assigned crustal values for these parameters; below this depth mantle values are 
assigned. For most of the models discussed here, we have chosen T c = t c . We do not mean to 
imply that thicknesses of the Martian lithosphere and crust coincide; this choice serves only to 
simplify the models. We have performed additional calculations with T e > t,, to explore the effect 
of a strong upper mantle on the flexural solutions. We take T c to be variously 20, 40, and 60 km. 
These values fall within the range obtained by Comer et al. [1985] for the Tharsis Montes and 
Olympus Mons. The outer radial boundary of the grid is taken as r max = 1200 km. 

As a check on the use of TECTON for plate flexure calculations, we have compared the 
stresses and displacements calculated by the code to those from an analytic solution for thick plate 
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flexure in axisymmetric geometry [Comer , 1983], We constructed a model of the plate only, 
without the volcano elements, similar to analytic models. This was done to avoid local stiffening 
of the plate due to the added thickness of the volcano near r = 0, which would change the resulting 
flexural profile. The load of the volcano was simulated by applying appropriate axisymmetric 
forces to the nodes on the top surface of the plate. The surface elevation of the plate at the final 
time step, relative to its value at the right boundary r = r^^ is shown in Figure 7a; also plotted is 
the analytic thick-plate solution. The difference between the two solutions is plotted in Figure 7b. 
The TECTON solution closely matches the analytic solution except for errors of less than 1% near 
r =0 and about 1.5% near the flexural arch. 

We note that for cases in which the volcano is solidly attached to the lithosphere, the deflection 
given by TECTON at the final time step is less than that of the analytic solution near r = 0. Radial 
distances of the zero crossing and flexural bulge are also greater than those for the analytic 
solution. These effects result from a stiffening of the plate by the volcanic load itself. Near r = 0, 
the elements that make up the volcano take up some of the horizontal compressive flexural stresses 
that would otherwise be taken up by the top half of the plate. This effect may be important when 
analytic thin-plate flexure models are compared to topography in order to estimate elastic plate 
thi cknesses, in that for welded volcanoes, T e would be overestimated. 

An important limitation of analytic models of lithospheric loading is that the load is emplaced 
instantaneously. This assumption implies a growth time t g much less than the mantle Maxwell time 
x M . For Mars, if T| ~ 10 21 Pa-s andp~ lO* 1 Pa, x M - 10^ s or 300 yr. A more realistic model 
would allow the volcano to grow incrementally with time, simulating the way in which many 
discrete lava flows are superposed to form an edifice. We have implemented a modification to the 
finite element code that allows us to specify the properties of a given element as a function of time. 
The finite element grid is set up exactly as in Figure 6a, except that volcano elements are initially 
given a density of zero and a low but finite value of Young's modulus (the latter to prevent a 
singularity in the stiffness matrix used to calculate displacements). These elements are considered 
"off;" that is, for the purposes of the simulation, such elements are considered initially unoccupied 
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by volcano material. At the beginning of an appropriate time step, an element can be turned "on" 
by changing the values of the' density and Young's modulus to normal crustal values. Here we 
have simulated the incremental growth of a volcano by switching on elements in the volcanic 
edifice in three stages (Figure 6b), at time zero, at 1000 Maxwell times, and at 2000 Maxwell 
times. We note here that our attempt to model a growing load is necessarily a coarse one. A 
typical terrestrial volcanic shield will consist of thousands of units emplaced over about 10^ to 10^ 
years (this timescale may be 10 8 to 10 9 years on Mars). Such units are much lower in volume than 
those we model. Further, our increments are separated by only 1000 Maxwell times (about 30,000 
years for the parameters adopted here). This value was selected to be long enough to complete the 
flexural response to a given increment, abut is not meant to imply that units are emplaced that 
rapidly. 

In order to explore the hypothesis that volcanoes on Earth and Mars may have similar tectonic 
structures, we have used special features in the TECTON code to simulate the effect of a 
detachment between volcano and lithosphere as described by Lipman et al. [1985] and others. The 
effects of fault slip along such a detachment on the stresses and deformation can be modeled using 
the slippery node method implemented in TECTON by Melosh and Williams [1989]. We apply 
slippery nodes to all the nodes originally on the line z = 0, i.e., along the interface between the 
volcano and the lithosphere. This modification allows the volcano to thrust outward over the 
flexed lithosphere. 

The evolution of the stress field within and beneath the volcanic edifice in a model is depicted 
by means of symbols for the principal stresses within each element (Figure 8, symbols after 
Melosh and Williams [1989]). The first and second symbols denote situations where the azimuthal 
or hoop stress is the intermediate principal stress &2- The third and fourth symbols represent cases 
where the hoop stress is the greatest extensional stress 03. The fifth and sixth symbols represent 
cases where the hoop stress is the greatest compressive stress Oj. 
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Given stress values, we use the Mohr-Coulomb failure criterion to estimate regions where 
faulting would occur. The Mohr-Coulomb failure equation relates shear stress x at failure to 

normal stress a n 

^failure = c + G n tan <j) (1) 

where c is the cohesive strength of the rock and + is the angle of internal friction. We adopt values 
for c (3.8xl0 7 Pa) and <> (49°) appropriate for basalt [Handin, 1966]. Elements where the shear 
stress has exceeded the Mohr-Coulomb criterion are shaded to indicate failure. Once we have 
found where failure is expected, the orientations of principal stresses are used to determine the 
style and orientation of faulting, according to the criteria of Anderson [1951]. Given the principal 
stress directions, the type of faulting (normal, thrust, or strike-slip) and the orientation (radial or 
circumferential) can be determined. Figure 8 gives types and orientations of faulting in elements 
with the given orientations of the stress symbols. This classification scheme holds if the principal 
stresses approximately correspond to the horizontal (a„ t Oqq ) and vertical stresses (c^ listed 
here. If the stresses are not aligned in this way, fault geometry may be less simply described. 
Strike-slip faults do not quite fit the usual definitions of radial and circumferential, because such 
faults are expected generally to strike obliquely to the principal stress directions. For most 
situations, strike-slip faults make smaller angles with the Ci direction than with the o 3 direction. 
Thus, when Gj = age we will consider the resulting faults approximately circumferentially 
oriented, and when Oj = c rr we consider them approximately radial. Care must be taken when 
applying these results to regions at depth, or where principal stresses are oblique to the surface. 

For a general state of stress in a plane layer constrained laterally and subjected to self-loading 
(before the volcanic load is applied), the orientation of the stress symbol will be as in the first line 
of Figure 8: maximum compression vertical and maximum extension horizontal. The magnitude 
of the deviatoric stress will increase with depth. This result can be derived from the equations for 
uniaxial strain (compaction) in the z-direction [e.g., Turcotte and Schubert , 1982 , p. 108]: 
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°rr=°ee=(f^)<*zz ( 2 ) 

For Poisson's ratio v = 0.25, the factor v/(l - v) is 1/3. Thus, the vertical stress will be three times 
as great as the horizontal stress, and with vertical stress equal to the overburden pressure (a^ = 
p c g Z , where p c is the density of crustal material and g is the gravitational acceleration), the stress 
difference will increase with depth. The state of stress in the terrestrial crust is still a matter of 
debate. McGarr [1989] argues that a stress state in which all principal stresses are equal to the 
overburden load is a better reference stress state than the "Poisson" stress state described by 
equation (2). However, it seems likely that the state of stress ties somewhere between the two 
extremes. In general, the stress state described by equation (2) is the initial state in the lithosphere 
for most of our models. We consider additional models, however, to determine the effects of an 
isotropic lithospheric prestress (with G n = Gqq = g zl = p c gz). For all models, the prestress in the 
asthenospheric mantle is taken to be isotropic. 

When applying the results of these calculations to actual geologic and tectonic features seen on 
Martian volcanoes, one must keep in mind the model limitations. The simplest calculations that we 
perform here start with the instantaneous placement, at time t = 0, of a significant portion of a 
volcanic edifice, equivalent to the condition t g « Ty . This assumption is valid only for 
individual flows constituting at most a minor mass fraction of the volcano. While we have 
modified the finite element code to accommodate changes in element properties with time, as noted 
above the individual load increments are large fractions of the volcano mass. Other limitations are 
also noteworthy. It is difficult to determine the geometry of a given volcano load, even at the end 
of the last major shield-building eruptions, because of subsequent flexure and deformation that 
may have significantly modified the volcano's characteristics. These calculations are performed 
under the assumption of axial symmetry; possible effects of nonaxisymmetric loading (such as 
regional stress) and complex three-dimensional geometry are not incorporated. The effects of 
magma pressure, transport, and evacuation are not addressed in these models; these processes 
likely have important influences for caldera and flank tectonic evolution [Thomas et al., 1990; 
Zuber and Mouginis-Mark, 1992]. Further, lateral variations in material properties (due to 
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horizontal temperature gradients, for example) are not included. Our choice that T e = t<., while 
made for convenience, is not necessarily valid for the Tharsis region, finally, it must be 
remembered that a viscoelastic model cannot account for relief of stress due to faulting (except for 
faults modeled a priori such as in the detached volcano models). Once faulting occurs, the stress 
fields calculated are no longer strictly valid, since faulting would relieve stresses locally and could 
thus alter the predictions for failure at subsequent times. As discussed below, some of these 
limitations can be relaxed in future modeling efforts. 

NUMERICAL RESULTS: STRESS, DEFORMATION, AND FAULTING 

Models with Growing Loads 

An illustrative model for the evolution of volcano-related stresses is shown in Figure 9, which 
depicts the case of an incrementally grown volcano on a lithosphere with T e = 40 km. At t = 0 
(Figure 9a) deviatoric stresses display the orientations and magnitudes expected from simple self- 
loading of horizontal layers everywhere but along the top layer of volcano elements, where the 
stresses are rotated such that the compressive axis is almost horizontal. As the effects of flexure 
manifest themselves with time, this rotation propagates deeper into the volcano, eventually 
reaching the upper part of the underlying plate (Figures 9b-f). These figures also show predicted 
locations and types of faulting as functions of time. For the initial volcano of 120 km radius, 
faulting is initiated in two areas by about 20 Maxwell times (Figure 9b). A subsurface region of 
normal faulting at a depth of 35-40 km (at the base of the elastic lithosphere) appears beneath the 
center of the volcanic load. Faulting also starts at the plate surface just outward of the volcano. In 
this region, all elements but one are predicted to experience circumferentially oriented normal 
faulting (the innermost element is predicted to fail by strike-slip). By 100 Maxwell times (Figure 
9c), both regions have expanded radially and vertically, and a new zone of faulting appears along 
the upper volcano flanks. The predicted fault style and orientation on the flanks is predominantly 
circumferential thrusting. Mostly radial normal faults are predicted for the subsurface failure 
region. By 1000 Maxwell times (Figure 9d), the area of circumferential normal faulting extends up 
to 300 km radius (not shown in figure) and up to 10 km depth. The area of strike-slip faulting has 
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grown in radial extent at the expense of the region of circumferential normal faulting, but each of 
the elements in this zone has already been predicted to fail by a circumferential normal mechanism 
at an earlier time step (see Figures 9b, 9c). The upper flanks of the volcano are predicted to fail by 
circumferential thrust faulting, and the lowermost element on the flank by radial thrust faulting. 

The failure region at the bottom of the plate now extends slightly beyond the radius of the volcano 
itself. Near r = 0, this region includes the entire bottom half of the plate (20 to 40 km below the 

original surface). 

In Figures 9e and 9f we see the effect of the next increment of volcano growth. Notice that 
most elements on the volcano flank predicted to fail in Figure 9d are no longer predicted to be in 
failure in Figure 9e, but some of them are predicted to fail again by the time depicted in Figure 9f. 
No failure is predicted, however, in the layer of elements constituting the new surface of the 
volcano. In Figure 9g, the effects of the final load increment can be seen. Only two elements in 
the entire edifice have deviatoric stresses sufficiently great to cause failure. 

The magnitude of the maximum deviatoric stress (the difference between the principal stresses 
of greatest and least magnitude) in the vicinity of the volcano is shown versus time m Figure 10. A 
comparison with the stress orientation symbols of Figure 9 indicates the principal stresses used to 
calculate the maximum deviatoric stress at a given point. Figure 10a shows the state of deviatonc 
stress of the elastic solution. The stress difference increases with depth, as discussed above. At a 
given depth, the deviatoric stress is higher beneath the initial volcano load increment than outward 
of the volcano. In Figure 10b (at t = 1000 x M ) concentrations of high deviatoric stresses have 
developed in the three characteristic regions of failure from Figure 9. Deviatoric stresses exceed 
150 MPa near the volcano summit and in a wide area near the bottom of the plate beneath the 
volcano. Failure is predicted in both areas. Deviatoric stresses are not quite so high in the near- 
surface region immediately surrounding the volcano, but faulting is predicted there because of the 
low confining pressure. There are two important regions of relatively low deviatoric stress; 
directly beneath the volcano summit near z = -4 to -14 km, and far from the volcano (at both the 
top and bottom of the plate). In Figure 10c (at t = 2000 t M ) the region of high stress in the edifice 
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has intensified, but remains centered at about the same position as before, at the top of the first load 
increment. The second load increment, emplaced atop the first, shows a much lower level of 
deviatoric stress than evident below. The deviatoric stress low beneath the volcano has shrunk in 
size and lowered its mean level to about z = -10 km. By the time subsidence is nearly complete 
(Figure lOd), most of the region near the volcano experiences deviatoric stress in excess of 200 
MPa. The regions of high deviatoric stress have expanded to maximum size, and the deviatoric 
stress minimum below the volcano center has again shrunk slightly. 

Changing the plate thickness can affect the evolution of regions of failure. For T e = tc = 20 
km , the entire flank is predicted to fail by thrust faulting after 1000 x M (Figure 1 la). 

Circumferential thrust faults on the upper flanks will be surrounded by a small annulus of radially 
oriented thrusts on the lower flanks. The region of failure outward of the volcano is smaller in 
radial extent than for the case with T e = 40 km (Fig. 10b). While strike-slip faulting in the inner 
portion of this region, initial predicted failure for most of these elements occurs by circumferential 
normal faulting. In any case, this inner annulus will be covered by volcanic material later in the 
growth of the volcano. At the final stage (Figure lib, r = 3000 x M ), the subsurface region of 
failure reaches through almost the entire thickness of the plate to the base of the volcano. At the 
radially distal edge of this region, stresses are rotated such that strike-slip faulting is expected, and 
the region merges with the failure region outside the volcano to form a continuous zone predicted 
to be faulting. More elements (5 vs. 2) in the edifice retain stresses large enough to cause failure 

than for the case with T e = 40 km. 

The evolution of maximum deviatoric stress for this case (Figure 12) follows a course similar 
to that for T e = 40 km. The regions of high deviatoric stress on the volcano flanks and at the 
bottom of the plate closely surround a stress minimum whose center lies at about z = -3 km. After 
the first load increment is accommodated (Figure 12a), peak stresses in the edifice are nearly twice 
as large as for the case with T e = 40 km. A comparison of the late-stage stresses for these two 
cases (Figures lOd and 12b) shows that the deviatoric stresses near the base of the plate are similar 
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in magnitude, although the vertical gradient of stress is greater for the thinner plate. Peak 
deviatoric stresses in the volcanic edifice are significantly higher for the thinner plate. 

Results for a model with a thicker elastic lithosphere and crust (T e = tc = 60 km) also reveal 
interesting changes from the first model. After the first loading increment (t — 1000 Fig* 13s), 
faulting has initiated at the lithosphere base and outside the volcano, but not in the edifice itself. In 
the outer faulting region we again observe elements where predicted failure evolves from initial 
normal to later strike-slip. The area of predicted strike-slip faulting is much larger in extent than 
for the cases with thinner lithospheres. After the final load increment (Figure 13b) the region of 
predicted near-surface normal faulting expands outward to 590 km from the center (Table 3), as 
well as down to a maximum depth of 20 km. The region of failure beneath the volcano is found at 
a deeper level (30-60 km depth), because of the greater plate thickness. Faulting of any type is 
absent from the edifice throughout the period modeled. 

Maximum deviatoric stresses for this case are shown in Figure 14. Compared to the case with 
T e = 40 km (Figure 10), there are several differences. At the final time steps shown (Figures lOd 
and 14b), the region of high stress near the base of the lithosphere (as defined by, say, the 300 
MPa contour) is greater in vertical and lateral extent for the case with T c = 60 km, even though the 
maximum stress deviator is about the same. The deviatoric stress minimum, beneath the volcano, 
is also larger in extent and centered more deeply (about -14 km) than its counterpart for T e — 40 
km. Deviatoric stresses within the volcanic edifice are significantly lower for this case compared 
with the other two. 

To investigate the effects of different crustal and lithospheric thicknesses, we performed a 
calculation with a 40-km-thick crust and a 60-km-thick elastic lithosphere. The upper mantle 
lithosphere has a higher value of Young's modulus and should cause the plate to behave more 
stiffly than if the plate consisted entirely of crustal material. Figure 15 shows plots of stress 
orientation and failure for this case; compare with Figure 9 for a case with crust and lithosphere 
both 40 km thick and with Figure 13 for a case with crust and lithosphere both 60 km thick. The 
evolution of the stress field and zones of predicted failure is qualitatively similar to those discussed 
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above. One effect of the stiff upper-mantle lithospheric layer is to concentrate into a narrower 
depth interval the zone of failure beneath the volcano, seemingly shielding the less rigid crust 
above from failure. In Figure 13, this zone extends from 30 to 60 km depth, whereas in Figure 15 
only elements in the mantle portion of the lithosphere (40 to 60 km depth) show stresses that 
satisfy the failure criterion. We also observe that the size of the zone of near-surface failure 
outward of the volcano, and the amount of strike-slip faulting within it, are both greater for the 
plate with a strong upper mantle lithospheric layer (Table 3). Contours of maximum deviatonc 
stress after the final load increment for this model (Fig 16) are broadly similar to Figure 14, except 
that the deviatoric stress level is much higher at the base of the lithosphere, revealing a 
concentration of stress in the strong mantle lithosphere. Also, the deviatoric stress minimum 

beneath the load occurs at a slightly shallower depth. 

Models with Instantaneous Loads 

It is instructive to consider a model with an instantaneously emplaced load, both as a 
comparison with analytic models and as a reference to study the effects that incremental volcano 
growth has on the evolution of stresses in volcanoes. Figure 17 shows stress orientations and 
failure for this case at the conclusion of flexure (t = 1000 x M ). These are qualitatively similar to 
those for the case with incremental loading (Figure 9g) except for the volcano flanks, which are 
predicted to have failed along the entire surface and to 15 km depth near the summit. 
Circumferential thrust faults are predicted on the upper flanks, and an annulus of radially oriented 
thrusts is predicted on the lower flanks. This state is similar to that predicted for the first load 
increment of a growing volcano (Figures 9d, 1 la). The maximum deviatoric stresses for this case 
(Figure 18) resemble those for the incrementally loaded case (Figure lOd) except for the 
concentration of stresses at the top of the edifice. In this respect the model behaves like an analytic 
plate flexure model, in which deviatoric stresses are a maximum at the top and bottom of the plate, 
with a minimum at mid-plate. The volcano serves effectively to stiffen the plate near r = 0, thus 
elevating the mid-plane from z = -20 km (for a 40-km-thick plate) to z = -8 (for a plate effectively 

65 km thick at r = 0). 
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Models with Detached Loads 

The evolution of stress orientations for an incrementally grown volcano detached from the 
lithosphere is shown in Figure 19. Immediately after loading, stresses in almost the entire edifice 
are oriented with maximum compression vertical and maximum extension out of the plane (Figure 
19a). This geometry may be compared with that in Figure 9a, where both the maximum and 
minimu m principal stresses are in the r-z plane in the upper flanks of the volcano, and the 
maximum compressional direction has rotated toward the horizontal. In Figure 19a the onset of 
faulting is predicted in two areas; at the plate surface just outside the volcano (as observed in most 
of the other models) and near the base of the volcano near r = 0. Radial normal faulting is 
predicted for this second region; the entire lower part of the volcano (near the detachment) exhibits 
horizontal extensional stress and vertical compressive stress, in contrast to the environment of 
generally horizontal compression in the models with welded volcanoes. This environment persists 
as the calculation proceeds (Figure 19 b-c). In Figure 19b, after 1000 Maxwell times, the area of 
radial normal faulting reaches outward through the edifice, above the detachment, to the surface at 
the volcano edge. Also note that stresses near the summit of the volcano have rotated to a state of 
horizontal compression. By the time the flexural response is complete (Figure 19c, 3000 Maxwell 
times), an annulus (with inner radius 170 km and outer radius 200 km) of radial normal faulting 
has formed on the lower flanks. Elsewhere, the stress state depicted in Figure 19c is similar to the 
final result of the welded case (Figure 9g), except that the area of predicted faulting outside the 
volcano is deeper, and the faulted area beneath the volcano does not extend as far vertically. These 
effects may be related to local stiffening of the plate by the welded volcano, which would raise the 

midplane of the plate. 

A plot of maximum deviatoric stress for this case (Figure 20) shows how the distribution of 
stresses is changed by the detachment. After first load increment, deviatoric stresses in the volcano 
tend to increase from top to bottom, the reverse of the pattern for the welded volcano (Figure 10b). 
A region of high deviatoric stress forms beneath the detachment (with approximately the same 
magnitude as the high stress region at the top of the edifice in Figure 10b). Figure 20b shows 
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stresses after completion of loading; stresses in the edifice are greatest at the bottom, directly above 
the detachment. The upper parts of the edifice have low levels of stress, except for a small area of 
slightly elevated stresses near the summit. The stress maximum beneath the detachment has grown 

in magnitude and area. 

In order to address the role of the initial stress state of the lithosphere in determining the time- 
dependent evolution of stresses and displacements, we consider a lithosphere with an isotropic 
prestress equal to the overburden pressure. For a model otherwise similar to that of Figure 9, the 
first faulting in the region outside the volcano (Figure 21a) is predicted to occur by strike-slip, in 
contrast to the prediction of circumferential normal faulting for the model of Figure 9. Notice that 
the onset of faulting in this region occurs at a later time, and the onset of faulting in the region 
below the volcano is delayed until about 400 x M (compared with 20 x M for the model in Figure 9). 
These delays are to be expected, since deviatoric stresses must build up from zero, as opposed to 
cases which start with a sizeable deviatoric stress. The extent of faulting at the final time step 
(Figure 21c) is similar (to 470 km radius, Table 3), and so is the predicted geometry of faulting in 
the region beneath the volcano (radial normal). In toe early stages of flexure, however, strike-slip 
faulting is predicted for a greater proportion of the near-surface region than for the model of Figure 
9 (Table 3). After the flexure following the final load increment the proportional area of predicted 
strike-slip to circumferential normal faulting at the surface is still about 10% greater than that of the 
model of Figure 9 (Table 3). Thus we conclude that the initial stress state of the lithosphere prior 
to volcanism can affect the distribution of near-surface faulting surrounding the volcano. 

Maximum deviatoric stresses (Figure 23) far from r = 0 are small, as might be expected in a plate 
which starts out with isotropic prestress. However, the stresses near the base of the plate beneath 
the volcano and in the volcanic edifice are qualitatively similar to those for cases with Poisson 

prestress. 

DISCUSSION: IMPLICATIONS FOR VOLCANO EVOLUTION 
The time-dependent displacement and stress fields calculated in our models can be used to 
predict the type and orientation of tectonic features that would result from loading of the lithosphere 
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by a volcanic edifice. We may compare these predictions to observed tectonic features on and 
around the Tharsis volcanoes. The evolution of tectonic deformation is also constrained by the 
model results, an advantage of the finite element method over analytical flexure and other 
lithospheric loading models which consider an instantaneously emplaced load. We may then 
evaluate die validity of each model and die applicability of terrestrial analogues for some tectonic 

and volcanic features. 

flank Tectonics: Circumferential Graben and Predicted Strike-Slip Faulting 
As noted earlier, circumferendal graben are observed on die lower flanks and the surrounding 
plains for all three of die Tharsis Montes (Figures 2-5). Formation of there circumferendal graben 
requires an environment of radial horizontal extension a, the top of the crust surrounding the 
volcano and in die near-surface region of die volcano flank. Almost even- mode, predicts 
circumferential norma, faulting in an annulus around the volcano. Models in which tire volcano .s 
welded to the underlying plate give near-surface horizontal compression in die edifice after 
lithospheric flexure, while models which include a detachment between volcano and plate show 
horizontal extension in the lower volcano flanks. A, die base of the detached volcano, die hoop 
stress Gee is extensional, and die radial normal stress <s„ is compressional but small compared with 
rhe vertical nonnal stress o. In order for circumferential extensions! features to form, 0[r should 
be the leas, compressive stress. However, since die magnitudes of a„ and a e9 are small in this 
region, a small perturbation due to local or regional nonaxisymmetric stresses may determine die 
orientation of failure features [Nakamura, 1977], The Tharsis Montes are located slightly below 
.be crest of die Tharsis Rise, on a broad slope to die northwest [U. S. Geoiagica, Survey 1989], 
This long-wavelength slope of Tharsis, in die presence of a detachment, would be expected to add 
a northwest-southeast extensional stress; superposition with the axisymmetric model stress field 
would result in radial extensions! features in the northeastern and southwestern quadrants (and, at 
sufficiently high stiains, a NE-SW bending rift), and circumferential graben on the northwestern 

and southeastern flanks. 
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On Arsia Mods, graben are prominent on the half of the volcano centered on the northwestern 
azimuth (all the way up the flank to the caldera rim. Figures 2, 5), but are generally absent on the 
upper southeastera-centered half. Several graben appear at the southeastern edge of the edifice, 
however. These graben are likely to be related to extensional stresses at the surface of the 
lithosphere surrounding the volcano, as described above, and thus not related to stresses in the 
edifice per se. The graben on the northwestern half terminate near the northeast-southwest- 
trending rift that bisects the volcano. We suggest that a portion of the northwestern half of the 
volcano was able to slide down the slope of the Tharsis rise, creating a stress state that led to the 
formation of the numerous graben. In the southeastern half, however, buttressed by the rising 
slope in that direction and thus unable to move, the stress state remained stable (horizontal 
compression at the surface, as in Figure 9f). This scenario would account for the dramatic 
difference in faulting character on the two "halves” of the volcano. Of course, stresses resulting 
from caldera formation and magma chamber foimation may have had a strong effect on the stress 

field near the summit as well [Zuber and Mouginis-Mark 1992]. 

We can apply the above arguments to the remaining Tharsis Montes. Graben on the flanks of 
Pavonis Mons are in general narrower than those on the northwestern flank of Arsia Mons, but 
more widely distributed in azimuth relative to Arsia Mons [Zimbleman and Edgett, 1991]. The 
widest graben are on the west and northwest flanks, analogous to the situation on the flanks of 
Arsia Mons. The graben on the north, northeast, and east flanks of Pavonis Mons, however, 
i-m-af that the stress state in the edifice is too complex to be described solely by the above 
scenario for Arsia Mons. Ascraeus Mons shows the least amount of normal faulting on its flanks; 
tire few concentrations of graben he on the lower slopes, far from the summit caldera. Thus 
Ascraeus Mons shows the least evidence of tectonic modification due to the effects of a 
detachment Perhaps the detachment is smaller in extent there because of inherent differences in 
crustal properties, or because it has had less time to develop. The latter scenario is consistent with 
the conclusion of Cnmpler and Aubele [1978] that Ascraeus Mons is at an earlier stage in its 
development than the other two volcanoes. 
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Stresses from both analytic flexure models and our models with instantaneous volcano loads 
predict strike-slip faulting at the surface of the plate, immediately outside the load Evidence for 
strike-slip faulting around any of the Tharsis volcanoes, however, is lacking. The absence of such 
features has been discussed by Golombek [1985] and Schultz and Zuber [1992], in the general 
context of large axisymmetric loads on planetary lithospheres. Examination of the plots of stress 
orientation and failure from our models (Figures 1 1, 13, 15) however, reveals that the predicted 
failure modes for elements in the region surrounding the volcano can change with time. In this 
region, every element that eventually fails (except the innermost) is predicted to do so first by 
concentric normal faulting. This can be seen more clearly in a Mohr diagram of the stresses in a 
typical element in this area (Figure 22). The Mohr circles grow in size as flexure proceeds, until . 
the failure envelope is exceeded at about 70 tm after first loading. The largest circle shows stresses 
at the conclusion of flexure. This stress state corresponds to that given by analytic plate models (in 
which flexure occurs instantaneously). Applying a failure criterion simply to the final stress state 
can be misleading. Consideration of the stress history in our models indicates that at the time of 
first failure, the shear failure criterion is satisfied in the surface elements immediately outside the 
volcano, and the principal stress directions predict circumferential normal faulting. These first- 
formed faults would tend to relieve extensional stresses in their vicinity and would serve as 
preferred planes of weakness for the relief of stresses accumulated later. Strike-slip features would 
thus tend to be prevented from forming. We also note that for the lower range of values of 
lithosphere thickness considered here (20-40 km), the zone of strike-slip faulting predicted from 
the late-stage stress models lies at the edge of the volcanic construct, which has the lowest surface 
elevation (Figure 7). This topographic low could accumulate lava flows and erosional deposits 
which would tend to cover deformational features (circumferential graben as well as strike-slip 
features). 

Horizontal Compression in the Edifice 

Circumferential terraces high on the flanks of Olympus have been interpreted as concentric 
thrust faults by Thomas et al. [1990], who applied elastic finite element models to investigate 
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edifice stresses and deformation. Our simplest model (a welded, instantaneously emplaced 
volcano) predicts a similar stress state, with circumferential thrust faulting on the upper flank and 
radial thrusts on the lower flank (Figure 17). Plate flexure models also predict large horizontal 
compressive stresses in the area of the load, which increase with height above the mid-plane of die 
plate, reaching a maximum at the surface, where thrust faulting is predicted. Models with growing 
loads, however, display a different distribution of stress. The maximum horizontal stress occurs 
in the lowest (first) load increment and decreases upward with each later load increment (Figures 9- 
14). Each successive load increment "feels" tectonic stresses from itself and all later load 
increments. The las. units emplaced are sensitive only to die effects of their own loading and 
therefore exhibit a low degree of horizontal compression. Thus, thrust faults (predicted by simpler 
models) associated with such stresses would be expected to occur either early in the evolution of 
the volcano or not at all and, if formed, would be covered by later units which remain unfaulted 
because *e failure criterion is not satisfied (Figures 9 . 1 1). We note, however, that stresses in 
these elements are close to failure. Local stress permrbafions and/or variations in internal fricuon * 
could result in some failure visible at die surface. On the other hand, if loading increments are 
smaller than those adopted here (which is likely), then dm degree of horizontal compression from 
such elements would be smaller than those shown here, and faulting would not be expected. 

The presence of graben high on the flanks of the Tharsis Montes is not consistent with the 
state of horizontal compression in die edifice found in the welded volcano models. Odier factors, 
such as a basal detachment, lithospheric rifting, or caldera stresses must be invoked to explain die 

presence of these faults on the upper flanks. 
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The volcanic edifice stress fields calculated here have important implications for magma 
transport and eruption. Magma propagates through the lithosphere along fractures that form 
perpendicular to the direction of leas, compressive stress. Immediately following emplacement, die 
stresses in a given load increment are characterized by o 3 horizontal (e.g.. Figure 9a). Thus show 
that nugma can propagate vertically to the summit region. As a result of flexure, however, the 
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principal stress directions in that load increment rotate such that (Ji is horizontal and (73 is vertical. 
Magma propagation should tend to be horizontal in such at layer, leading outward from the summit 
to the flanks. Stresses in earlier load increments remain in horizontal compression as volcano 
growth proceeds. This behavior results in a maximum in horizontal stress near the core of the 
volcano, with lessening horizontal compression at higher (later) load increments. This stress 
pattern contrasts with that from analytic flexure calculations, where stresses increase linearly with 
distance from the plate midplane. Thus an implication of the models treated here is that within an 
interval of shield building there should be an evolution in favored eruption location from summit to 
lower flank. Such an evolution is consistent with the first two stages in the sequence of major 
events for the Tharsis Montes determined by Cnanpler and Aubele [1978], 

The stress evolution of a growing volcano may affect the location of magma collection zones 
or chambers. The maximum horizontal compression is achieved within the earliest load 
increments, near the bottom of the load. Such a region of high compressive stress at lower levels 
in the edifice could inhibit ascending magma from reaching the summit, and could thus lead to the 
pooling of magma. Since magma will propagate perpendicular to the direction of least compressive 
stress, radial propagation of magma in sheet dikes or sills might also occur. On the other hand, 
high stresses in this region will be relieved by faulting, until the failure criterion is no longer 
satisfied. This region may thus be the most pervasively faulted and fractured zone, and thus 
perhaps a favored location for magma to collect. 

The issue of compressional versus extensional horizontal stresses in the volcanic edifice is 
critically important when considering volcano growth. In order for the volcano to grow, there 
must be a path for magma to reach the surface (or near the surface in the case of growth by 
intrusion). Large horizontal compressive stresses within the edifice may cut off summit eruptions 
unless there is sufficient magma overpressure. Thus, the route from a central magma chamber to 
the summit may be cut off, and flank eruptions would be favored, as noted above. After sufficient 
flank eruption, however, loading on the flanks might then lead to relief of compressive stresses 
near the summit, re-opening the path to the summit caldera. When the effects of time-dependent 
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loading during volcano growth are included, the magnitude of horizontal stresses in the later stages 
a of evolution are greatly reduced compared with stresses at the equivalent level in instantaneously 

pvAH models When both time-dependent loading and a basal detachment are included, horizontal 
; compressive stresses in the edifice are reduced ftrrther, and a state of extensional hoop stress is 

reached in the lower part of the edifice, adjacent to the detachment (Figure 19). Summit eruptions 

5 are then more easily accommodated. The stress state in the lower edifice resembles that determined 

seismically for western Hawaii by Gillard et al [1992], with a, vertical (or nearly so), and both <r 2 
- mi oj horizontal and similar in magnitude. The model stress state also resembles that of a tectonic 

model for the base of Hawaii proposed in Figure 3 of Wyss and Koyanagi [1992] (based in part on 
the results of Gillard et al [1992]). The model of Wyss and Koyanagil 1992] has <r 3 parallel to 
die direction of slip along the basal detachment, whereas our model has 03 perpendicular to the 
basa) slip. However, since <72 and 03 are similar in magnitude [ Gillard et al. 1992], local 
perturbations and departures from axisymmetty may determine the orientations of these stresses in 

actual volcanoes. 

Flank Tectonics: Aureoles, Lobes and Scarps 

Z Results of our models also help us to understand the flank tectonics of theTharsis volcanoes. 

in particular the large downslope motions of material in the fonn of large slides or slumps. As 
— noted above, the Olympus Mons aureole has been interpreted as a large landslide deposit resulting 

from slip along large-offset listric normal faults [Francis and Wadge. 1983], the upper reaches of 
which are postulated to coincide with the basal scarp. Aureole material is most prominent to the 
north and west; there is an aureole lobe to the southeast, however, suggesting that the Tharsis rise 
did not fully buttress the southeastern flank of Olympus. Evidence for similar mass movement 
activity on a smaller scale is seen on the northwestern flanks of Ascraeus and Arina Montes 
[Zimbleman and Edgett, 1991]. 

As discussed above, the interaction between the axisymmetric stress field from the volcano 
loading and regional stresses due to topography may influence the location, orientation and extent 
of scarps and of landslides triggered by fault sUp on the Tharsis volcanoes. The Olympus Mons 

1 
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aureole deposits are of greatest extent to the north and west, but also there is a significant aureole 
lobe to the southeast. The basal scarp has greatest relief on the northwestern and southeastern 
flanks. Models with a detachment predict that at the edge of the edifice CT 3 = Oqq and is 
extensional. This can account for radial normal faulting but not for the scarps and aureole 
deposits, which in the above scenario require that 03 = G„. The addition of a regional northwest- 
southeast extensional stress, induced by the regional gradient in elevation, however, would 
provided the necessary orientation of stresses in the northwestern and southeastern quadrants of 
the volcano, where the aureole and scarp are most prominent. The aureole deposits to the north 
and west of Olympus Mons could then be a result of large landslides analogous to those off the 
Hawaiian chain but on a larger scale. The basal scarp would then consist of the coalesced head 
scarps of these landslides. The lobate features to the northwest of Arsia and Pavonis Montes 
would then be similar in origin but on a smaller scale, without prominent basal scarp formation. 
This scenario also explains why the deposits are always more prominent on the northwestern side 
of these volcanoes, since this is the regional downslope direction. Note that the above scenario 
postulates a consistent stress field (of NW-SE extension) throughout the evolution of the volcano. 
This contrasts with the view of Francis and Wadge [1983], who suggest that proto-Olympus Mons 
was elongate to the northwest due to an early NE-SW extensional regime, which allowed 
preferential effusion along a NW-SE axis. 

Large landslides and slumps, triggered by extension on the flanks, may be a second important 
mode of volcano growth, in that mass movement clears the way for more volcanic material to be 
emplaced behind it [Upnum et al . , 1985]. The model of Lipman and others emphasizes the role of 
riilfft intrusion at the rift zones, with the subvolcano detachment passively taking up the resulting 
strain. Our results suggest that without such a detachment, the rift zones may not be able to form, 
due to a predicted environment of horizontal compression. As a welded volcano grows, an ever 
increasing horizontal compressive stress near the base of the edifice would eventually choke off 
magma supply to the summit. Thus, the Hawaiian rift zones may owe their existence to the 
presence of the detachment, and volcanoes lacking such a detachment may exhibit important 
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structural and evolutionary differences from Hawaiian volcanoes. For example, Nakamura [1980] 
has suggested that the lack of rift zones in the Galapagos shields is attributable to a paucity of 
oceanic sediments on the young ocean floor beneath the volcanoes. 

Linear Rift Zones 

Linear symmetries observed in all the Tharsis Montes volcanoes may be indirect evidence for 
failure in the sub- volcano region. The second stage in the Crumpler and Aubele [1978] sequence 
of volcano evolution involves the development of a linear rift bisecting the volcanoes, with 
eruptions emanating from this rift. Radially oriented normal faulting under a volcano might be 
thought at first to tend to divide the volcano into radial sections, much like pie slices. The 
existence of the linear rifts on the Tharsis Montes, however, suggests that this mode of failure may 
be concentrated into one large linear feature that bisects the volcano. This feature may owe its 
orientation to regional stress fields or to pre-existing zones of weakness. Our calculations predict 
radially oriented normal faulting in a broad zone beneath a volcano. For the case with T e = 20 km 
(Figure 1 1), this zone extends beyond the volcano radius and reaches quite close to the surface. 
Such a zone, if concentrated into a linear feature and if continued to the surface by sustained 
faulting, may give rise to the linear rifts and bilateral symmetry observed on the Tharsis Montes. 
Stresses due to topography of Tharsis might be expected to add a component of northwest- 
southeast extension (the direction of greatest slope). This would concentrate the radial faulting into 
a rift zone striking northeast-southwest, as is observed on the Tharsis Montes. The rift zones 
would form orthogonal to the regional topographic gradient, as is observed. A problem with this 
model is that the rifting has to propagate upward somehow to the surface of the lower volcano 
flanks. The prediction of radially oriented normal faulting on the lower flanks for the model with a 
sub-volcanic detachment (Figure 20) suggests another mechanism for forming the observed rifts. 
This explanation has the advantage that rift formation would occur early, almost immediately after 
volcano creation, a timing more consistent with the chronology of Crumpler and Aubele [1978]. If 
the detachment model is applicable to Mars, it may help explain the existence of the rifts. 
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The evolution of Pacific hotspot volcanoes may provide some insight into the formation and 
influence of a bisecting rift zone. For instance, Holcomb [1985] has suggested that the island of 
Molokai in the Hawaiian chain resembles a volcano cut almost exactly in half, with one part 
transported downslope in one or more massive submarine slides, and the other part remaining to 
form the present island. Bathymetry of the ocean floor north of Molokai [Moore, 1964] supports 
this conclusion. We have earlier conjectured that a similar but less extreme mechanism could 
explain the presence of graben on the upper slopes of Arsia Mons and Pavonis Mons, in that as the 
volcano halves slide apart from each other, an extensional environment is created near the summit. 
To investigate this possibility quantitatively, a fully three dimensional model, or an axisymmetric 
model with non-axisymmetric loading, is required. 

Detachment Tectonics 

The flank tectonic characteristics of the Tharsis volcanoes, combined with the observations of 
and modeling results for terrestrial volcanoes with detachments, strongly suggests that detachments 
analogous to those of Hawaii-type volcanoes [Lipman et al., 1985] exist on Mars. Such a 
detachment cannot arise from the presence of a thick oceanic sediment layer as for Hawaii. It has 
nonetheless been demonstrated that liquid water existed in moderate quantities at the Martian 
surface earlier in the history of the planet [e.g., Carr , 1987]. While liquid water is now unstable 
at the surface, large quantities of water are believed to be buried as ground-ice. Tanaka [1985] 
proposed that gravity sliding and spreading could account for the Olympus Mons aureole (and 
basal scaip) providing that a basal detachment zone of material with low shear strength was 
present. He calculated that a layer of material with about 10% interstitial or interbedded ice could 
provide the necessary detachment. The pre-volcamc upper crustal layer may have consisted 
dominantly of an impact breccia regolith produced by repeated impacts. Such a heavily fractured 
material is likely to be porous and thus may contain interstitial ice. Even in the absence of an 
excess ice content, the mechanical weakness of such a layer could contribute to the formation of a 
detachment between a strong lower crust and a volcanic edifice. Tanaka [1985] states that such a 
structure is unique to Olympus Mons. However, Timbleman and Edgett [1991] have mapped 
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lobate deposits on the northwestern flanks of all three Tharsis Montes that may be smaller-scale 
analogues to the Olympus Mons aureoles, so some form of detachment beneath all four constructs 
appears to be indicated. 

Plate Thickness and Volcano Size 

The ratio of elastic plate thickness to volcano diameter is observed to have significant effects 
on the evolution of the subsidence, stress, and regions of predicted failure in our computations. 

We note several trends as plate thickness increases relative to volcano diameter. Study of the 
patterns of stress orientations and predicted failure in Figures 9, 1 1, and 13, shows that the outer 
region of failure, at a given time, increases in horizontal (and vertical) extent, and a greater 
proportion of the inner part of this region is predicted to fail at later stages by strike-slip faulting, 
with increasing plate thickness (Table 3). The deviatoric stress plots of Figures 10, 12, and 14 
show that the average depth and volume of the region of low deviatoric stress beneath the volcano 
also increase with increasing plate thickness. The magnitude of deviatoric stress within the 
volcanic edifice, and in particular the magnitude of maximum horizontal compressive stress, 
decreases with increasing plate thickness. As noted above, a zone of large horizontal compressive 
stress may affect the propagation of magma from its mantle sources to the volcanic edifice. Elastic 
plate thickness may thus play a role in determining the location of magma storage areas. For the 
largest plate thickness considered here, stresses in the edifice never exceed the failure limit (Figure 
13). 

The relative timing of events in our models is generally consistent with the chronology set 
forth by Crumpler and Aubele [1978] for the Tharsis Montes. After construction of the main 
shield, rift formation and flank eruptions along this rift occur in the second stage of evolution. 
Figure 9b shows that, after 10 Maxwell times, the principal model result is the beginning of 
rotation of the stress axes in the edifice (which is expected to favor flank eruptions). After a few 
tens of Maxwell times the failure zone beneath the volcano (which could favor formation of the rift) 
has begun to grow. The thir d stage involves summit subsidence and formation of concentric 
fractures and graben. After 100 Maxwell times, our computations show increased subsidence as 
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well as the growth of the region of normal faulting outward of the volcano (Figure 9c-g). The 
fourth stage of Crumpler and Aubele [1978] involves further eruption along the rift, flooding the 
summit and flanks. Late stage flank eruptions are consistent with the stress fields in our 
calculations. Summit eruptions, however, would not be favored by the computed stress fields. It 
is likely that the stress field is influenced by the effects of previous stages of the evolution. Flank 
eruptions result in an added flank load, not modelled here, to which the lithosphere will respond 
with additional flexure. Subsidence and additional flexure near the flanks of the volcano may at 
least partially relieve the compressional stresses near the summit, thus enhancing the possibility of 
additional summit eruptions. By this scenario, the preferred site for eruptions could alternate 
between summit and flank, on time scales on the order of 100-1000 Also, the effect of 
sections of the volcano sliding away from the bisecting rift may result in an extensional 
enviro nm ent near the summit, once again favoring summit eruptions. 

Future Directions 

Because of the limitations of these calculations in accurately modeling the behavior of large 
volcanoes, we hope to be able to relax some of these limitations in future models. In our models, 
assume that the crust and mantle are perfectly viscoelastic; within an element, no provision is made 
for stresses that exceed reasonable failure limits (such as in all the shaded elements in Figures 9, 
11, and 13). The inclusion of localized yielding in a volcano model can be presumed to have 
important effects on the geometry and magnitudes of the calculated stress and deformation fields. 
A finite element code that incorporates plastic yielding will allow us to account more realistically 
for the effects that failure will have on predicted fault locations and styles. Further, the effect of 
thermal stresses has been ignored. Thermal stresses can play an important role in the evolution of 
the volcano, both through stresses accumulated by cooling after emplacement of each layer, and 
through the time-dependent thermal stresses in regions reheated by the later passage or storage of 


magma. 
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CONCLUSIONS 
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A mass movement origin has been suggested for the vast aureole deposits surrounding 
Olympus Mons [Harris, 1977; Lopes et al. , 1980; Francis and Wadge, 1983] as well as smaller 
lobate deposits to the northwest of each of the Tharsis Montes [Zimbleman and Edgett, 1991]. 
Analogy with large slumps and slides off of the Hawaiian islands [Moore et al., 1989; Lipman et 
al., 1985] suggests that a detachment between volcano and lithosphere may play a role in the 
evolution of volcanoes on Mars as well as Earth. Models that include such a detachment, in 
contrast with models with welded volcanoes, predict a state of horizontal extension in the edifice 
(and eventually radial normal faulting on the lower flanks). When the effects of non-axisymmetric 
stresses due to regional topography are considered, predicted patterns of faulting match well 
observed structures on the Tharsis Montes. 
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FIGURE CAPTIONS 

Figure 1. Viking Orbiter view of Olympus Mons. The multiple-pit caldera and basal scarp (at 
bottom) are the main tectonic features visible. No graben are seen on the flanks of the volcano 
or on the volcanic plains immediately surrounding the structure. Frame V0641A52; width of 
image is 440 km. 

Figure 2. (a) Viking Orbiter view of western Ascraeus Mons. The complex caldera is seen at the 
s ummi t, and graben with a predominantly circumferential orientation are seen along the 
western and southwestern margins. Graben are located on Ascraeus Mons itself and in 
surrounding flows (unit as. Figure 3); some graben grade into collapse pits or sinuous 
channels formed by magma withdrawal and downslope flow. Frame V0643A78; width of 
image is 420 km. (b) Viking Orbiter mosaic of the western flank of Arsia Mons [Eliason et al. 
1991]. Note that graben extend all the way up the flank to the edge of the summit caldera, (c) 
Viking Orbiter mosaic of the eastern flank of Arsia Mons [Eliason et al. 199 1]. Note the lack 
of graben and the radial texture of flows on the upper flank. Graben are evident near the 
boundary between the edifice and the surrounding plain. 

Figure 3. Geologic map of Ascraeus Mons and surroundings, simplified from Scott et al. [1981b] 
by Comer et al. [1985]. Volcanic units shown include relatively young Ascraeus Mons flows 
(as), intermediate-age Tharsis Montes flows (tm), and volcanic material undivided by age (vu). 
Also shown as a distinct unit is slide material (s), interpreted by Scott and coworkers as 
landslides and debris flows. Dashed lines show approximate elevation contours, in kilometers, 
relative to a fourth-degree, fourth-order equipotential [Wu, 1978]. The summit caldera 
complex is indicated by inward hatched lines. Extensional faults and graben are shown as 
heavy lines. 

Figure 4. Geologic map of Pavonis Mons and surroundings, simplified from Scott et al. [1981a, 
b, c] and Scott and Tanaka [1981] by Comer et al. [1985]. Units shown, in addition to those 
described in Figure 3, are relatively young volcanic flows from Pavonis Mons (pm) and Arsia 
Mons (am). Other information follows the format of Figure 3. Circumferential graben reach 
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quite far up the slopes, and an approximate bilateral symmetry can be seen about an axis 
trending approximately NNE-SSW. 

Figure 5. Geologic map of Arsia Mons and surroundings, simplified form Scott et al. [1981c] by 
Comer et al. [ 1985]. See Figures 3 and 4 for further explanation of symbols. Circumferential 
graben extend almost up to the summit caldera complex. 

Figure 6. (a) The finite element grid used for the calculations shown in later figures. The volcano 
is the small triangular area in the upper left-hand comer, (b) Plot of element types for models 
with incremental volcano loading, implemented by means of time-dependent element 
properties. Elements marked with a "3" are "on" from the start of the calculation; those with a 
"4" are switched "on" after 1000 Maxwell times, and those with a "5" are switched "on" after 
2000 Maxwell times. 

Figure 7. Comparisons of the TECTON solution and analytic thick-plate flexure [Comer, 1983] 
for vertical displacements resulting from an axisymmetric load simulating a conical volcano. 
For both cases, T e = tc = 40 km. (a) surface elevation relative to that at r = r^ after 
essentially complete viscoelastic relaxation in the asthenosphere; (b) difference between the two 
solutions. 

Figure 8. Key for stress symbols in models. These symbols denote the principal stress directions 
under axisymmetry. The hourglass shapes are oriented along the direction of greatest 
compressive stress, the bars along the direction of least compressive stress. Magma is 
expected to ascend locally along the plane orthogonal to the least compressive stress, that is, 
perpendicular to the bar. 

Figure 9. Principal stress orientations and predicted failure in the vicinity of an incrementally 
grown volcano load on a 40-km-thick elastic lithosphere. See Figure 8 for the meaning of the 
stress orien tati on symbols. Elements in which stresses satisfy the Mohr-Coulomb failure 
criterion are shaded. Given the orientation of the stress symbol in a shaded element, the style 
and orientation of faulting can be inferred, (a) The elastic solution, (b-d) 20, 100, and 1000 
Maxwell times after the first load increment is added; (e-f) 20 and 1000 Maxwell times after the 
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second load increment is added; and (g) 1000 Maxwell times after the third load increment is 
added. 

Figure 10. Contours of maximum deviatoric stress in the vicinity of the volcano, for the model of 
Figure 9. Contour interval is 50 MPa. (a) The elastic solution, (b-d) 1000 Maxwell times 
after the first, second, and third load increments are added, respectively. 

Figure 11. Principal stress orientations and predicted failure for an incrementally grown volcano 
load on a lithosphere with T c = 20 km. (a) 1000 Maxwell times after the first load increment is 

added; (b) 1000 Maxwell times after the third load increment is added. 

Figure 12. Contours of maximum deviatoric stress for the model of Figure 1 1. (a) 1000 Maxwell 
times after the first load increment is added; (b) 1000 Maxwell times after the third load 
increment is added. 

Figure 13. Principal stress orientations and predicted failure for an incrementally grown volcano 
load on a lithosphere with T e = 60 km, 1000 MaxweU times after the third load increment is 

added. 

Figure 14. Contours of maximum deviatoric stress for the model of Figure 13, 1000 Maxwell 
times after the third load increment is added. 

Figure 15. Principal stress orientations and predicted failure for an incrementally grown volcano 
load on a crest 40 km thick and on an elastic lithosphere 60 km thick, 1000 Maxwell times after 
the third load increment is added. 

Figure 16. Contours of maximum deviatoric stress for the model of Figure 15, 1000 Maxwell 
times after the third load increment is added. 

Figure 17. Principal stress orientations and predicted failure for an instantaneous volcano load on 
a lithosphere with T e = 40 km, after 1000 Maxwell times. 

Figure 18. Contours of maximum deviatoric stress for the model of Figure 17, after 1000 
Maxwell times. 

Figure 19. Principal stress orientations and predicted failure for a model with a growing volcano 
load, a basal detachment between lithosphere and volcano, and crest and elastic lithosphere 40 
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km thick, (a-b) 10 and 1000 Maxwell times after the first load increment is added; (c) 1000 
Maxwell times after the third load increment is added. 

Figure 20. Contours of maximum deviatoric stress for the model of Figure 19. (a) 1000 Maxwell 
times after the first load increment is added; (b) 1000 Maxwell times after the third load 
increment is added. 

Figure 21. Principal stress orientations and predicted failure for a growing volcano load with 
isotropic prestress in the lithosphere. Crust and elastic lithosphere are 40 km thick, (a-b) 40 
and 1000 Maxwell times after the first load increment is added; (c) 1000 Maxwell times after 
the third load increment is added. 

Figure 22. Mohr diagrams for typical elements in regions of predicted failure. Diagonal line is the 
Mohr-Coulomb failure envelope for basalt [Handin, 1966] with an angle of internal friction <J> = 
49° and cohesion c = 38 MPa. (a) Mohr diagram for an element at the top of the lithosphere 
(stresses calculated at element center r — 265 km, z = -2.5 km). In order of increasing radius, 
the circles represent the stress state at t = 0 (elastic solution), 20, 100, 1000, 2000, and 3000 
tm • When the failure envelope is first exceeded (r = 70 x m ), 03 = Orr and a i * a zz- At the 
final time shown, a 3 = a„ and 0\ = 099 . See Figure 9 for the principal stresses used to 


calculate circles. 
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TABLE 1. Dimensions of Tharsis Volcanoes 


Volcano 

Diameter, 

km 

Relief, 

km 

Olympus Mons 

600* 

24» 

Ascraeus Mons 

400 b 

18 b 

Pavonis Mons 

320 b 

14 b 

Arsia Mons 

420 b 

19 c 


*WuetaI. (1984). 
b Blasius and Cutts (1976). 
c Blasius and Cutts (1981). 
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TABLE 2. / 
Parameter 

Adopted rare 
Crust 

lUiClVi * 

Mantle 

E, Pa 

lxlO 11 

3xl0 n 

p, kg/m 3 

3000 

3500 

V 

0.25 

0.25 

Lithosphere 

Asthenosphere 


t\. Pa s 


1x1027 1x1021 
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G = growing volcano, 

W = welded volcano, 

1 = instantaneously emplaced volcano, 

D = detached volcano, 

PS = isotropic pre-stress applied to lithosphere, integer = elastic lithosphere 
thickness. 

2 Radius values given at I =■ 1000 tm after application of given load increment 

3 Value less than the maximum radius of the volcano at a given load increment 
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Arsis 2.2 13 40km 1.66700BH1 seconds 

stress orientation/Failure plot full output timestep # 4 

cohesion = 3.50000Bf07 pa; phi - 49.0000degrees 
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Arsia 2.2 13 40km time= 8.33400EU1 seconds 

stress orientation/failure plot ^ output timestep # 8 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Arsia 2.2 13 40km time* 8.33300B+-12 seconds 

stress orientation/failure plot frU out P ut timestep # 12 

cohesion* 3.50000Bf07 pa; phi * 49.0000degree$ 
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Figure 9d 

























































































































































































Arsia 2.2 13 40km time= 8.50000E+12 seconds 

stress orientation/failure plot frU output timestep # 14 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Arsia 2.2 13 40km — 

stress orientation/failure plot full output timestep # 

cohesion = 3.50000E+07 pa; phi = 49.0000degrees 
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Figure 11a 
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Figure 13a 
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Figure 15 












































































































































































































time = 8.33300E+12 seconds 

Arsia 2.2 11 40km 

. w full output timestep# 12 

stress orientation/failure plot 

cohesion = 3.50000B+07 pa; phi = 49.0000degrees 
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Arsia 2.2 13 slip 40km «— 8.333005,10 seconds 

stress orientation/failure plot full output ttaestep# 2 

cohesion = 3.50000Ef07 pa; phi = 49.0000degrees 






Figure 19a 
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Figure 19c 
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Features on Venus Generated by Plate Boundary 
Processes 
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successfully used to understand the tectonics of the Earth can also be 
applied to Venus. 


tuc r arrro - ECLOGUE PHASE TRANSITION AND THE 

Technology, Cambridge, MA 02139. 
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Among the four mountain belts surrounding w . t oWhmi Planum. The 

similar to those of terrestrial tholeutic basalt [10]. Because the rp\o \\ the 7-1 1 km 

OTsfb teen 

thermal boundary layer is assumed to occur y P^^ s ^ thermal diffusivity is assumed 

fixed at 750 K and at a value T« taten “ *“ fSfite basalt [13], and the densities of gabbro 

■ 5* 



essure increases at a kayw* 
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well as phase changes. The volume fraction of reacted component, y, is given y, 

\jnj/ = D/r 2 

D = Dai Opx = LlxlO 5 exp (-400 kJ/RT ) 

wW R the eas constant and T is the absolute temperature. The diffusion rate, however, is 

«perimentaUy uncertain^xrause - 

lithosphere, v is obtained by integration over time. The density at a given ucp. 
“SS' Temperances KddST^s. and manfle .cMri ferrates of 


ELEVATION OF MOUNTAIN BELTS ON VENUS: Namiki N. and Solomon S.C. 


horizontal convergence of 1 <H 5 (Fig. la) and 1<W* s' (Fig. lb). ForjU 
fTable U thicknesses of crust and thermal boundary layer are assumed to be initially 2 U ana ou 
£L respectively and to increase to 100 and 250 km, respectively. Temperature profiles for the 
Sain rateof Iff 1 ’ 5 s 1 are vertically stretched as the crust and lithospheric mantle are thickened 
(Fig la). Temperatures do not increase significantly from initial values because heat 
transferred by advection and the contribution of crustal heat production is minor. 
remains meltable for 50 My or more, and the elevation of shoyt^ htihos^wcM ino^» as 
much as 12 km above the surrounding plains in that time intend (Fig. 2a). ^c phax umsition 

proceeds, i.e., elevation is limited, only if grains are small crust 

p For *j=io-i6 s -i t ^ 5^1 heat production dominates advective heat transfer after me crest 

becomes as thick as 60-80 km (Fig. lb). The resulting increase in temperature hastens die phase 

SS SI Si also lengthens rite formation time of 

characteristic reaction time, which depends on r. D, and temperature at the base of crest. For 
brew t^ns(r=f0 mm), the elevation reaches U km or more if the tnthal tempera^ at 

of tite crust. T is 1 150 K (Model 6 in Fig. 2b). For the same value of T c the elevation is at most 6 
tata^tinlofliun raZ (Model 1). This result constitutes ^Wer bound on T e for small 
if n_r)_ r , j s assumed, that upper bound is lowered to 1050 K (Model 8). 
^Discussion. Because at long wavelengths the topography of western Ish^- Terra 
with the gravity field dynamical support of the broad 4-km elevation of the region is tuceiy ie.g., 

S of Maxwell Montes above the adjacent plateau is a more 

^a^SlSns^TSn maximum relief than tile 1 1 -km elevation above mean planetary radius. 

The results are insensitive to the assumed initial thickness of crest but are sensitive to die density 
w^n cras t and mantle. If densities of 3000 and 3300 kg nr’ for crust and mant e are 
asrerredfekvations are -40% lower than the values presented here. Such changes ^nstratn the 
^Sof themial struaure and phase transition depth more severely at low strain rate than the 

de ”Two diffusion rates have been assumed in this study so as to represent ith 

diffusion data. We should note, however, that under wet conditions, i.e., in the presence of eimer 
water [16] or melt [17], grain-boundary diffusion becomes much more efficient vduree 
HiSncion This is potentially noteworthy for understanding the contrast between Maxwell and 
rj Tutnntes Despite the fact that Danu Montes display compressional deformation as extensive 
STeotet^^ elevation Fs i lit tie" as 0 km above tire tandq jtei. 

Such comparatively modest elevation may be related to the presence of magmatic features within 
Danu Montes if elevation is limited by an enhanced diffusion rate because of the melt at gran 
boundaries. Assessing the cause of higher temperatures beneath 
detailed thermal models than the simple one-dimensional model considered here. 

Conclusions Taking into account the temperature-dependent reaction rate of the gabbro 
edtSSSSSmdK horizontal strain rates of 10-« and 10“ resun in sigmffcant djfacoeei 
b th?of maximum elevation of mountains, not only because of rhti iifference in ^ 

for relief, but ^^^^^.i^SS^mti^.mountata 

belts can be explained as a consequence of disequilibrium phase teurvi^ depth te a mdejange 

of ohvsical parameters, although a comparatively young age for Maxwell Montes (50 y) 

F^liS FofSsrer horizontal strain rate of 10 -‘ s >, only linuted parameter values for thermal 

models are allowed. 
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Figure 1. Theraial evolution of crust 
and mantle thickened by a uniform 
horizontal strain rate of (a) 10- 15 and 
(b) 10* 16 s* 1 . The phase diagram of 
the gabbro (G) - gamet granuUte 
(GG) - eclogite (E) transition [13] is 
also shown. The models shown 
correspond to models (a) 1, 2 and 3, 
and (b) 4, 5 and 8 in Table 1. 
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for the models in 
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Abstract 

Maxwell Montes, standing up to 7 km above the adjacent highland plateaus, constitute the 
highest mountain belt on Venus. Because the thickness of the oust is likely to be limited by the 
gabbro - garnet granulite - eclogite phase transitions, this relief is difficult to reconcile with the 
assumption of thermodynamic equilibrium and a standard Airy isostatic model. We explore the 
hypothesis that the crust-mantle boundary is not in phase equilibrium, but rather is rate-limited by 
the temperature-dependent diffusion of the slowest ionic species. Under the simplifying 
assumption that the mountains formed by uniform horizontal shortening of the crust and 
lithospheric mantle at a constant rare, we solve the one-dimensional thermal evolution problem. 

The time-dependent density structure and surface elevation are calculated by assuming a 
temperature-dependent reaction rate and local Airy isostatic compensation. For a rate of horizontal 
strain of ICC 15 s 1 or greater, the rise in temperature at the base of the crust during mountain 
fomrim is modest to negligible, the deepening lower crust is metastable, and surface elevation 
increases as the crust is thickened. For strain rates less than 10 16 s' 1 , in contrast, crustal 
temperature increases with time because of internal heat production, and the lower crust is more 
readily transfomted to the dense eclogite assemblage. For such models a maximum elevation is 
reached during crustal shortening. While this maximum relief is 7 km or more for some models, a 
smaller density contrast between crust and mantle than assumed here (500 kg m‘3) and 
incorporation of horizontal heat transport would lessen this value. We therefore favor formation of 
tire mountain belt at a strain rare at least of order 10-15 s -l. B y this reasoning, MaxweU Montes 
must be comparatively young, of order 50 My. 


' "N 
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Introduction 

The linear mountain belts of western Ishtar Terra on Venus are notable for their topographic 
relief and slope and for the intensity of surface deformation [Barsukov et al, 1986; Solomon et al., 
1991, 1992]. These four mountain belts, including Akna, Freyja, Maxwell, and Danu Montes 
(Figure 1), are generally regarded as products of large-scale compression and horizontal shortening 
of die crust and lithosphere [Crumpler et al., 1986; Solomon et al., 1991, 1992; Kaula et al, 

1992]. Among these mountain belts, MaxweU Montes is the highest and stands up to 1 1 km above 
mean planetary radius and 7 km above the highland plain Ukshmi Planum (Figure 1). Volcanism 
is rare to absent in MaxweU, Freyja, and Akna Montes, but a number of magmatic features are 
evident in Danu Montes [Solomon et al., 1991; Head et al., 1991; Kaula et al., 1992], the 

mountain range least elevated above Lakshmi Planum. 

The bulk composition and radioactive heat production of the crust on Venus, where measured, 
are similar to those of terrestrial basalts [Surkov et al., 1984, 1987]. Because basalt transforms to 
eclogite at high pressure and because eclogite is likely to be denser than mantle material, the 
thickness of the crust may be limited by the gabbro - garnet granulite - eclogite phase transitions 
[Anderson, 1981; Turcotte, 1989]. The phase transition depth depends on the temperature 
structure of the crust beneath the mountains. If the thermal gradient is low and thermodynamic 
equilibrium is assumed, the phase changes would take place at shaUow depth, and simple isostatic 
models underpredict the retief \V order Bruegge and Head, 1991]. If the thermal gradient is high, 
in contrast, the base of the crust is predicted to undergo melting before a sufficiently deep root can 
form to support the mountains. Thus the 7-1 1 km maximum elevation of MaxweU Montes is 
difficult to understand under the assumptions of local isostasy and thermodynamic equilibrium 
except in the unlikely situation that the crust beneath the mountains contains a large volume of 
magma [Vorder Bruegge and Head, 1991]. 

A possible explanation for the great relief of MaxweU Montes is that the base of the crust is not 
in phase equUibrium. Because of the high surface temperature on Venus (750 K) and the very low 
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water abundance of the lower atmosphere [von Zahn « at.. 1983], (he Venus crus, is (hough, (o 
contain negligible water [Koala, 1990). I( has been suggested lhai under anhydrous condinons, 
the gabbro - eclogite phase transition takes place by solid-state diffusion and may require a 
geologically significant time to run to completion [Ahrens and Schubert. 19751. Solid-state 
diffusion is a strongly temperature-dependent process, so that quantifying this suggestion for 
application to Venus requires the solution of a heat transport problem. In this paper we develop 
simple models for die thermal evolution of die crust beneath mountain beta on Venus in an attempt 
to constrain the time-dependent depth of the gabbro - eclogite transition and thus to assess this 

hypothesis. 

Because at long wavelengths the topography of western Ishtar Terra is correlated with the 
gravity field, at least partial dynamical support of the regional elevation of about 4 km is likely 
[e.g., Grimm and Phillips. 1991]. We therefore adopt the 7-kro elevation of Maxwell Montes 
above the adjacent plateau as a conservative estimate of the relief attributable to local isostasy. 


THERMAL MODEL 


We assume that the crustal and mantle portions of the thermal boundary layer shorten 
horizontally and thicken vertically with time in the manner of pure shear. The thermal structure is 
then governed by the one-dimensional heat equation, 

(1) 

9t dz dz 2 Cp 


where T is temperature, . is time, i is die horizontal strain rate, z is depth, K is the dtemtal 
diffusivity, A is the crustal heat production, and Cp is the specific heat Under the assumption that 
defonnation is by pure shear, f is constant and the vertical velocity is given by die product of t and 
z. We solve equation (1) numerically by an explicit finite difference approximation. 


Model Parameters 

The dtemtal diffusivity is assumed to be 1.0xl0-‘ m* s > in both crust and mantle. The crustal 
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heat production is assumed to equal l.OxlO* 10 W kg* 1 , on the basis of the K, U, and Th 
concentrations determined by gamma ray spectrometry by Venera and Vega landers [Surkov ex al ., 
1987]. The specific heat is assumed to be 850 kJ kg* 1 K* 1 . The other physical properties of 
crustal and mande materials are assumed to be those, respectively, of tholeiitic basalt and peridotite 

[Basaltic Volcanism Study Project, 1981]. Temperature at the surface and the bottom of the 

thermal boundary layer are fixed, respectively, at 750 K and at a value T w , a free parameter in 
numerical models. The initial temperature profile is determined by the assumption of steady-state 
conditions with zero strain rate. The increase in temperature with depth due to adiabatic 
compression is not considered here, because compressional heating is negligible compared with 

crustal heat production. 

Tbi is interpreted as the temperature in an isothermal core of a thermal convection cell [Turcotte 
and Schubert, 1982] and is, thus, fixed at a constant For a given Ty, an initial temperature T c at 
the base of the crust is determined. Since temperature at the base of the crust controls the elevation 
of the mountains, T c has more importance than Ty. Zuber [1987] analyzed wavelengths of 
tectonic features on Venus and found that the crustal thickness is less than about 15 km if the 
average vertical thermal gradient in the crust dT/dz is 25 K km* 1 , and less than about 30 km if 
dT/dz is 10 K km* 1 . The absence of significant viscous relaxation of impact crater relief also 
limits the crustal thickness in plains regions to be less than 10 km for dT/dz = 20 K km* 1 and 20 
km for dT/dz = 10 K km* 1 [Grimm and Solomon , 1988]. These results give a maximum 
temperature difference across the crust of about 400 K. Therefore Ty is given initially so that T c 
satisfies this upper bound. 

The phase diagram is assumed to be that of tholeiitic basalt [Ito and Kennedy, 1971], and the 
densities of gabbro and eclogite are taken to be 2900 and 3500kg m* 3 . The density of garnet 
granulite is assumed to increase linearly from that of gabbro to that of eclogite as pressure increases 
at a given temperature. The density of the mantle is assumed to be 3400 kg m* 3 . The elevations 
calculated by an assumption of local Airy isostatic compensation are sensitive to the density 
difference between crust and mantle. We have adopted a relatively large value of 500 kg m 3 . 



Because a large density difference results in a higher elevation for a given root thickness, this 
assumption permits a conservative lower bound on horizontal strain rate. We discuss later the 
consequence of a lesser density contrast 

Diffusion Rates 

The micromechanism of the gabbro - eclogite transition is not well understood. In this study 
we assume that volume diffusion of cations is the process most likely to limit the rate of the 
transformation, which involves chemical as well as phase changes. The volume fraction of reacted 
components, \| /, is given by, 

W = D/r 2 (2) 

where r is the typical grain radius and D is the diffusion coefficient [Ahrens and Schubert, 1975]. 
For each parcel of shortening thermal lithosphere, y is obtained by integration over time. The 
density at a given depth is determined from the volume fractions of reacted and unreacted 
components. 

There are several reactions occurring among constituent minerals in the garnet granulite stability 
field [Ringwood, 1975]. The rate of each reaction is difficult to assess, because of a paucity of 
experimental data. We focus on the diffusion rates of ionic species in pyroxene, garnet, and 
olivine, because those minerals characterize the phase assemblages and therefore play important 
roles in the phase transitions. Plagioclase is possibly another important mineral for this 
assessment, but no data for self diffusion in plagioclase are available. The data for volume 
diffusion in silicates show that silicon and oxygen diffuse much more slowly than divalent and 
trivalent cations. Thus the Si-O groups provide a static framework through which the cations 
diffuse [Freer, 1981]. We examine the diffusion rates of four major cations, Mg 2+ , Fe 2+ , Ca 2+ , 
and AP+ (Figure 2). The slowest diffusion rate is likely to limit die reaction rate. 

Diffusion rates of Fe 2+ and Mg 2+ in garnet and olivine [Chakraborty and Ganguly, 1991; 
Morioka and Nagasawa, 1991], of Mg 2+ in pyroxene [Cygan and Lasaga, 1985], and of Ca 2+ in 
olivine [Morioka and Nagasawa, 1991] have been experimentally determined (Figure 2). Of these 



7 


the diffusion rate of Mg* to game. (D M g.Gt) is toe lowest Of course, if toe diffusion rate of any 
unmeasured major species is slower than that of Mg* to garnet then toe phase change rate will be 
limited by this slower cation. We toereforo take toe diffusion coefficient D„ g ,o. f« this process as 
an upper bound on D in equation (2). 

Other diffusion rates, while not well-determined experimentally, have been estimated, e' 
from measurements at a single temperature, or from observed compositional gradients. While all 
such estimates include large uncertainties, we take the minimum value of these estimates as a 
bound on D in equation (2). This minimum is for diffusion of A1 in ortoopyroxene (Figure 2), 
estimated from an analysis of compositional gradients in naroral assemblages [Smith and Barron, 
1991). The estimated rate depends on assumptions of cooling rate and pressure. Whileneitoer 
dam nor estimates are available for toe diffusion coefficients of Al* in olivine and gamed cation 
diffusion to olivine and garnet is generally faster than to pyroxene [Freer, 1981; Snath and Brown, 
1988], and thus is not likely to be rate-limiting. Hie diffusion of Ca* to game, has been estimated 
by numerical simulation only at 1200°C [Chakraborty and Ganguly. 1991] and is faster than that of 
Al* in ortoopyroxene at that temperature (Figuro 2). The mobilities of Ca* and Fe* are held to 
be comparable in game, and pyroxene, because neither game, nor pyroxene in natural assemblages 
is consistently homogeneous in those elements [Smith and Barron, 1991]. These considerations 
reinforce toe view that Da,.opx gives a reasonable lower bound on D in equation (2). The tivo 

bounds on D are 

D £ Dm.ci * 2-8 x 10-* exp [- (270 kj + 3.2 x 10- 4 P )/( RT )] (3d 

DSDMOrx-llkl^expt-dOOkJART)] (3b) 

where P is in Pa [Morioko and Nagasawa, 1991; Smith and Barron, 1991]. 

There is some possibility that toe assumed lower bound on D is too large. Sautter oral. (1988] 
and sooner and Hone (1990] reported a diffirsion rate of Al* to clinopyroxene less than toe lower 
bound given above (Figure 2), while Freer et al. [1988] experimentally obtained a maximum 
diffusion rate as groat as toe upper bound given above (Figure 2). Furtoer, toe diffusion rates to 
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plagioclase are generally less than in pyroxene [Smith and Brown, 1988], and interdiffusion of 
CaAl-NaSi in plagioclase is several orders of magnitude slower than the lower bound on D [Grove 
et al., 1984], If D is as low as suggested by these last results, no phase transition is expected to 
occur during the formation of mountains on Venus. However, the reaction rate is unlikely to be 
much slower than as given by the lower bound on D, because if volume diffusion is too slow, the 
reaction rate will be limited by grain-boundary diffusion. 

Grain-boundary diffusion, in fact, has been proposed as a possible governing process in 
metamorphic reactions, because grain-boundary diffusion coefficients can be several orders of 
magnitude greater than those for volume diffusion [. Joesten , 1991]. Whether the phase change 
proceeds by volume diffusion or grain-boundary diffusion will be governed by the process that is 
most efficient for mass transport. Grain-boundary diffusion coefficients can be misleading, 
however, because mass transport by such a process is confined to a thin layer, i.e., the grain 
boundary. Unfortunately, grain-boundary diffusion data for the minerals involved in the gabbro - 
garnet granulite - eclogite reactions are few and generally not well determined [Joesten, 1991]. 
Reliable experimental data for grain-boundary diffusion of oxygen in forsteritc, for instance, have 
been reported only at two temperatures [Condit et al., 1985]. These measurements suggest that 
grain-boundary diffusion for grain sizes of 1-10 mm is no more efficient for bulk transport than the 
lower bound on volume diffusion adopted above (Figure 2). 

NUMERICAL RESULTS 

Temperatures in the thickening crust and mantle are calculated for rates of horizontal 
convergence of 10"^ (Figure 3a) and 10*^ s** (Figure 3b). For all models discussed here (Table 
1), thicknesses of crust and thermal boundary layer are assumed to be initially 20 and 50 km, 
respectively, and to increase to 100 and 250 km, respectively. Temperature profiles for the strain 
rate of 10- 15 s 1 are vertically stretched as the crust and mantle portion of thermal boundary layer 
are thickened (Figure 3a). Temperatures do not increase significantly from initial values because 
heat is transferred mainly by advection and the contribution of crustal heat production is minor. 


Hence gabbro remains metastable for 50 My or more, and the elevation of shortened lithosphere 
can increase as much as 12 km above the surrounding area of undeformed lithosphere in that time 
interval (Model 1 in Figure 4a). Equation (2) shows that the reaction runs to completion in a 
shorter time for small grains than for larger grains. Also it is clear from equation (3) that a higher 
temperature increases the diffusion rate and promotes the reaction. For the slower diffusion rate 
(Daj.giX however, a strain rate of lfr 15 s' 1 is sufficiently fast to suppress these effects (Model 2). 
The phase transition proceeds, i.e., elevation is limited, only if grains are small, diffusion is rapid, 
and temperature is high (Model 3). 

For a strain rate of 10* 16 s' 1 , crustal heat production dominates advective heat transfer after the 
crust becomes as thick as 60-80 km (Figure 3b). The resulting increase in temperature hastens the 
phase transition, and the slower strain rate lengthens the formation time of the mountains. For 
larger grains (r=10 mm) and a diffusion constant given by equation (3a) the elevation reaches 1 1 
km or more even if T c is as great as 1 120 K (Model 4 in Figure 4b). For the same value of T c the 
elevation is at most 7 km for grains of 1 mm radius (Model 5). If T c is greater than 1 120 K, the 
model gives a maximum elevation lower than 7 km, that is, the model fails to explain the observed 
relief of Maxwell Montes. Therefore this result provides an upper bound on T c for small grain 
size. If D is given by equation (3b), that upper bound is lowered to 1030 K (Model 6). Thus the 
combination of the longer formation time due to the lesser value of y and the faster reaction due to 
higher temperature yields constraints on the range of parameters controlling the phase change if the 
elevation of the mountains is to reach 7 km or more above the adjacent plateaus. 

DISCUSSION 

The thermal model described above is simple, but significant constraints on the strain rate 
governing mountain budding and the temperature at the base of the crust are nonetheless obtained. 
In this section, we first consider how the results would differ if key parameters were varied or the 
models were modified. We then discuss some implications of die numerical results given above. 

The thermal evolution results are insensitive to the assumed initial thickness of the crust. At a 
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strain rate of 10-15 r l, to advection-dominated thermal profile indicates that the temperature of 
each parcel of crust changes little as the parcel deepens. For a given value of T c , the thermal 
structure at a given crustal thickness is not affected by initial crustal thickness, or, eqmvalendy, 
Wtial thermal gradient Only the time required for the crust to reach that thickness dependsonthe 
initial thickness. At a strain rate of 10-“ tr*. in contrast the thermal proSe is close to steady state, 
that is, the geotherm is nearly the same as that for zero strain rate and thus docs not depend on the 
initial crustal thickness. The initial crustal thickness is still important for the maximum relief, 
however, because by assumption it is identical to the crustal thickness beneath surrounding areas, 
taken here to be undergoing negligible rates of horizontal strain. By the assumption of Airy 
isostasy, at a given crustal thickness beneath the mountains a difference of 10 km in the thickness 
of the crust beneath surrounding areas results in a 1.5 km difference tn relief. 

Also important is the density difference between crust and mantle. If densities of 3000 and 
3300 kg nr 3 for crust and mantle are assumed, relative elevations are -40% lower than the values 

prasentedhere. At a strain rate of lO-lSs-l.gabbroiclowercrustremainsmetastableduringthe 

formation of mountains in most cases, as shown by models 1 and 2. Therefore a smaller density 
contras. Simply requires thicker crush l.e„ greater duration of shortening, and does not otherwise 
constrain the model parameters. In contrast, a. a strain rate of Iff* r 1, maximum elevations are 
limited in most cases by the gabbro - eclogite phase transition because of tire high temperature at 
the base of the oust, as shown by models 5 and 6. The lower relief resulting form a smaller 
density contrast lessens the upper bound on T c . Thus such a lower density contrast constrains the 
models of themral structure and phase transition depth more severely at low strain rales than do tire 

density values adopted above. 

At a strain rate of 10->« r ' >, the thermal regime is dominated by crustal heat production and 
diffusion of heat Horizontal heat transfer within the shortening thermal lithosphere therefore is 
probably not negligible. The maximum increase in temperature at the bottom of the crast due fo 
horizontal heat transfer can be estimated form the solution for the instantaneous heating of a semi- 
infinite half-space [Turcotte and Schubert, 1982], 
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5T = AT eTfc (4) 
whOT 8T is the increase in temperature, AT is the temperature difference benealh die crust and an 
isothermal core which surrounds the shortening lithosphere, L is the horizontal scale of the 

, ..ins and t is the formation time. In our problem, AT is given by Tm - T c , and L is -100 

fan. Formation times of 50 My at a strain rate of 10-« s> and of 500 My at a strain rate of 1(H« 
s-l result in 5T of 2x10^ and 60 K, respectively. The assumption of one-dimensional heat transfer 
is therefore reasonable at the higher strain rate, but horizontal diffusion of heat increases 
temperature significandy at the lower strain rate. If horizontal heat transport were included in the 
thermal models at the lower strain rate, the reaction rate of the phase change would be faster, the 
maximum elevation would be less, and the corresponding upper bound on T c would be lower. 

These considerations favor the formation of Maxwell Montes at horizontal strain rates on the 
order of 10* 15 s* 1 . The implied age of the mountains is of order 50 My or less. Such an age is 
young compared with the average crater retention age of the surface of Venus of about 500 My 
[Phillips et al., 1991]. An important question is whether such a young age is consistent with the 
presence of the 100-km diameter Cleopatra crater [Basilevsky et al, , 1987] (Figure 1) in MaxweU 
Montes. The age of a single crater, of course, even a large one, is difficult to specify. The rate of 
formation of impact craters larger than 20 km in diameter on Venus has been estimated as 3.3 ± 1.8 
x 10- 15 km- 2 yr 1 from the statistics of Earth-crossing and Venus-crossing asteroids and an 
assumed set of scaling laws [Shoemaker et al ., 1991]. If MaxweU Montes, 1000 km long and 500 
km wide, formed in 50 My, the expected number of craters larger than 20 km in diameter is -0.08. 
The probability that one and only one crater larger than 20 km diameter formed in MaxweU Montes 
in 50 My is thus ~8 %. This value indicates that the presence of a large impact crater is not likely if 
the mountains formed in only 50 My, but the figure is not so low as to reject the hypothesis at high 

confidence. 

Although Danu Montes display congressional deformational features as extensive as the other 
mountain belts of Ishtar Terra, their maximum relief with respect to Lakshmi Planum is as Utile as 
1 km (Figure 1). The presence of magmatic features within Danu Montes [Solomon et al ., 1991, 
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1992; Head et al., 1991; Kaula et al., 1992] indicates that temperature beneath the mountains 
exceeded the solidus of crustal or shallow mantle material at least locally. Given such evidence for 
relatively high temperature in the lower lithosphere, the reaction of sufficiently deep lower crustal 
material has probably gone nearly to completion, which may account for the comparatively modest 
topographic relief of these mountains. It is also noteworthy that even a small amount of melting 
can greatly enhance grain-boundary diffusion [Condit et al. t 1985; Watson, 1991]. Assessing the 
cause of higher temperatures beneath Danu Montes requires more detailed thermal models than the 
simple one-dimensional ones considered here. , 

CONCLUSIONS 

Taking into account the temperature-dependent reaction rate of the gabbro - eclogite phase 
transition, horizontal strain rates of 10 15 and 10“ s 4 result in significant differences in the 
maximum elevation of mountains on Venus, not only because of the difference in the formation 
time for relief, but also because of the difference in the thermal regime from advection-dominated 
to crastal-heat-production dominated. For a strain rate of lO 15 s 4 , the observed maximum 
elevation of Maxwell Montes can be explained as a consequence of local isostasy and 
disequilibrium phase boundary depth for a wide range of physical parameters. For the lesser 
horizontal strain rate of 10 46 s 4 , only limited parameter values for thermal models are allowed, 
and this limitation becomes more difficult to satisfy when a smaller density contrast between the 
crust and mantle is assumed or when horizontal heat transport is considered. We favor the 
formation of Maxwell Montes at a strain rate of 10 45 s 4 or greater, although the likelihood that 
Cleopatra crater is younger than the implied age for Maxwell Montes (-50 My) is low. Magmatism 
within Danu Montes indicates higher temperatures in the crust or shallow mantle beneath those 
mountains than beneath the other mountains. The modest elevation of Danu Montes relative to the 
other mountain belts of Ishtar Terra can be understood as a consequence of near completion of the 


phase transition. 
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FIGURE CAPTIONS 


Figure 1 . Topographic map of western Ishtar Terra. ConretnmapofMagellan altimetry data 
t&foi a, 0.5° intervals of latitude and longitude; 1-tan contour interval; polar stereographic 
projection. Cross symbol shows the location of Cleopatra crater (65.9“N 7.0‘E). 


Figure 2 Diffusion constants for ionic species in olivine, pyroxene, and garnet Solid tines show 
diffusion constants that were experimentally determined (1 is for Fe in game. [Ckakraboriy and 
Ganguly, 1991]; 2 is for Mg in garnet [Chakraboriy and Ganguly, 1991]; 3 is forFe in oltvme 
[Morioka and Nagasava, 1991]; 4 is for Mg in olivine [Marioka andNagasava, 19911; 5 is for 
Mg in clinopyroxene [Cygan and iaraga, 1985 ]; and 6 is for Q in olivine [Morioka and 
Nagasawa, 1991]). The shaded area shows the range of estimates for Daia». including 
maximum and minimum C7) values [Smith and Barron, 1991). Cress symbols show the values 
determined a. a single temperamre (8 is for Cain game. [CHakraboriy and Gang*. 1991]; 9 ts an 
upper bound on the diffusion constant of A1 in clinopyroxene [Freer et of.. 1988]; 10 and 1 1 are 
values for A1 in clinopyroxene fiom Sautter et al. [1988] and Sander and Hone [1990], 
respectively). The constants for grain-boundary diffusion of O in foreterite [Condi, e, al., 1985] 
are converted re equivalent volume diffusion by multiplying by tite ratio of grain boundary width re 
grain diameter (12); the grain boundary width was taken re be 3 pm [Condi, e, al., 1985], and the 

grain diameters were taken to be 1-10 mm. 


figure 3. Thermal evolution of crus, (solid) and mantle (dashed) thickened a. a uniform horizontal 
strain rare of (a) 1(H ! and (b) l(h«s-i. The phase boundary of die gabbro (G), game, granulire 
(GG). and ec.o^re (E) assemblages and the solidus (S) (fro and Kennedy, 197 1] are shown by 
dotted lines. The models shown correspond to models (a) 1 and (b) 4 and 5 in Table 1. 


figure 4. Temporal variation of elevation for tite thermal models in Table 1. 
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TABLE 1. Parameters for Thermal Models 



Model 

r, nun 

Y. s’ 1 
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T b i,K 
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CRUSTAL DEFORMATION ON VENUS: Simons M. ct al. 


Results. In all our models, convection produces horizontal compressional stresses in 
lithosphere above downwelling mantle and extensional stresses in lithosphere above 
upwelling mantle. As the convective vigor increases so does the magnitude of the stress. In 
models with constant-viscosity mantle overlain by a constant-viscosity crust, stress in the 
crust reaches values in excess of 100 MPa in less than 100 My. We find that the rate of 
increase in compressive stress decreases with increasing crustal viscosity. This is because 
the stronger the crust, the more the development of the convective instability in the mantle 
driving the deformation is impeded. We also find that the magnitude of the peak compressive 
stress achieved above the downwelling increases with higher viscosities and/or with thinner 
initial crustal layers; the stronger the crustal lid, the more are tractions from mantle 
convection supported in the crust. Since force balance on the crust requires that shear 
traction integrated along the base be balanced by normal tractions integrated through its 
thickness, the thinner the crust, the larger the horizontal stresses. 

Both analytical models [9,10] and our numerical models of convection-induced crustal flow 
indicate that the amplitude and sign of the topography are highly time- and rheology- 
dependent. In general, possible responses of the crust to mantle flow can be divided into 
three categories. The first involves little, if any, crustal flow, and topography results mainly 
from the transmission of normal tractions induced by density contrasts within the mantle. 

The second possible regime involves substantial crustal flow, with geologically rapid 
thickening over convective downwelling and thinning over convective upwelling. In this 
regime the effects of crustal thickness variations dominate the topography. A third possible 
regime lies between the first two, with "in phase" deformation on short time scales and 
crustal flow on longer time scales. A strong mantle lithosphere tends to shield the crust from 
convective shear tractions, and topography results mainly from the transmission of convective 
normal tractions. A relatively weak lower crust facilitates crustal deformation, and the 
isostatic effects of crustal thickness variations dominate the topography. 

Consideration of geoid to topography ratios (GTRs) can restrict which regime of crustal 
response is appropriate for Venus. The distibution of estimated GTRs for several highland 
regions on Venus is bimodal with two clusters around 10 and 25 m/km [5]. The positive 
correlation of long-wavelength gravity and topography implies that there are no major regions 
that have negative GTRs. To keep the GTRs positive, the geoid must follow the surface 
topogaphy in sign at all times. Although the regime of negligible crustal flow can correctly 
predict both the sign and magnitude of the GTRs, it does not allow for the crustal deformation 
(i.e., flow) inferred from observations of tectonic features on Venus. 

In contrast to the regime of negligible flow, that of time-dependent crustal flow generates 
topography that changes sign. Over a mantle downwelling, the topography is negative in the 
early stages of deformation and positive in the later stages of deformation; the converse holds 
over a mantle upwelling. During the transitional period the topography goes through zero, the 
geoid does not go through zero, and the GTR is unbounded; this singularity is not observed 
on Venus. An example of the evolution of the GTR for a model in the time-dependent regime 
is shown in Figure 2. Note that before 90 My both topography and geoid are negative, but at 
90 My crustal thickening effects begin to dominate and the topography switches sign, causing 
the GTRs to first become unbounded and then negative. 

In the absence of a mechanism by which the sign of the geoid anomaly mimics that of the 
topography over a given upwelling or downwelling, only the regime of rapid crustal flow is 
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plausible. In the case of mantle down welling this would also require that the lower mantle be 
more viscous than the upper mantle in order to produce the required positive geoid anomalies. 
This has already been shown to be true for the Earth, where the observed geoid highs over 
regions of mantle upwelling and regions of mantle downwelling are best explained by the 
presence of a strong lower mantle [11,12]. The large positive GTRs and the presence of 
large shield volcanoes in certain highland regions on Venus, such as Beta Regio and Eistla 
Regio, are best explained as areas of mantle upwelling [5,13,14]. The regime of rapid crustal 
flow predicts crustal thinning over the upwelling. However, the extensive partial melt and 
ensuing volcanism expected over such regions of mantle may outweigh the effects of crustal 
thinning on the surface topography and thus also yield postive GTRs [15]. 



Figure 2. Evolution of GTR with time. At each time step, the GTR shown is the slope of the 
best fit line to the collection of geoid and topography values upward continued to 250 km 
above the surface nodes in the model. The result is for a downwelling model, with a crust of 
constant viscosity 50 times that of the mantle. In this model the sign of the topography is 
time-dependent, but that of the geoid is not. The singularity at about 90 My corresponds to 
the time when topography crosses the zero datum. Note that the GTRs for this model are 
neither always positive nor always bounded. 
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195, 1990; [9] D. L. Bindschadler and E. M. Parmentier, J. Geophys. Res ,95, 21, 1990; [10] 
H. Schmeling and G. Marquart, Geophys. Res. Lett , 17, 2417, 1990; [1 1] M. A. Richards and 
B. H. Hager, J. Geophys. Res , 89, 5987, 1984; [12] B. H. Hager, J. Geophys. Res , 89, 6003, 
1984; [13] R. E. Grimm and R. J. Phillips, J. Geophys, Res., 96, 8305, 1991; [14] R. E. Grimm 
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Geoid, Topography, and Convection-Driven Crustal Deformation on Venus; 
Mark Simons*, Bradford H. Hager*, and Sean C Solomon 2 , *Department of Eruth, Atmosphenc, 
and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 0213V, 
^Department of Terrestrial Magnetism, Carnegie Institution of Washington, Washington, DC 
20015. 

Introduction. High-resolution Magellan images and altimetry of Venus reveal a wide range of 
styles and scales of surface deformation [1] that cannot readily be explained within the classical ^ 
terrestrial plate tectonic paradigm. The high correlation of long-wavelength topography and gravity 
and the large apparent depths of compensation suggest that Venus lades an upper-mantle low- 
viscosity zone [2-5]. A key difference between Earth and Venus may be the degree of coupling 
between the convecting mantle and the overlying lithosphere. Mantle flow should then have 
recognizable signatures in the relationships between the observed surface topography, crustal 
deformation, and the gravity field [6,7], Therefore, comparison of model results with 
observational data qan help to constrain such parameters as crustal and thermal boundary layer 
thicknesses as well as the character of mantle flow below different Venusian features. We explore 
in this paper the effects of this coupling by means of a finite element modelling technique. 

Model The crust and mantle in the models are treated as viscous fluids. We solve both the 

equations of motion and the heat equation at each time step using a modified version of the two- 
dimensional, Cartesian finite-element program ConMan [81. Our formulation assigns material 
properties to each element and explicitly tracks the free surface topography and the crust-mantle 
interface by letting the grid deform with time in a semi-Lagrangian fashion. Hus procedure 
permits variation in the crustal buoyancy as well as specific crustal and mantle rheologies that can 
depend on composition, temperature and stress. In addition to the flow field, the stress field in the 
lithosphere, the surface topography, and the resulting geoid are readily calculated. The model 
domain is a square box. We impose frec-slip boundary conditions on the bottom and side walls, 
with no flow through these walls. The top of the box is a true free surface, so there is no need to 
derive surface topography from vertical stresses on the top of the box. We investigate two classes 
of models. In the first class flow is dominated by concentrated upwelling (Figure la), and in the 
second class flow is dominated by concentrated downwelling figure lb> We vary the initial 
crustal and thermal boundary layer thicknesses as well as the effect of crustal and mantle viscosities 
that are either constant, temperature-dependent, or fully non-Newtonian. 

Results. In all of our models, convection produces horizontal compressional stresses in 
lithosphere above downwelling mantle and extensional stresses in lithosphere above upwelling 
mantle. As the convective vigor increases so does the magnitude of the stress. In models with 
constant-viscosity mantle overlain by a constant-viscosity crust, stress in the crust reaches values 
in excess of 100 MPa in less than 100 My. We find that the rate of increase in compressive stress 
decreases with increasing crustal viscosity. This is because the stronger the crust, the more the 
development of the convective instability in the mantle driving the deformation is impeded. We 
also find that the magnitude of the peak compressive stress achieved above the downwelling 
increases with higher viscosities and/or with thinner initial crustal layers; the stronger the crustal 
lid, the more are tractions from mantle convection supported in the crust Since force balance on 
the’ crust requires that shear tractions integrated along the base be balanced by normal tractions 
integrated through its thickness, the thinner the crust, the larger the horizontal stresses. 

Both analytical models [9,10] and our numerical models of convection-induced crustal flow 
indicate that the amplitude and sign of the topography are time- and rheology-dependent In 
general, possible responses of the crust to mantle flow can be divided into three categories. The 
first involves little, if any, crustal flow, and topography results mainly from the transmission of 
normal tractions induced by density contrasts within the mantle. The second possible regime 
involves substantial crustal flow, with geologically rapid thickening over convective downwelling 
and thinning over convective upwelling. In this regime the effects of crustal thickness variations 
dominate the topography. A third possible regime lies between the first two, with "in phase" 
deformation on short time scales and crustal flow on longer time scales. A strong mantle 
lithosphere tends to shield the crust from convective shear tractions, and topography results mainly 


CRU! 


STTAL DEFORMATION ON VENUS: Simons M. et al. 


SSlSSSiS SSSSS! 3SS .SSSSffiS. *■ 

isbimodal with two clusters around 10 and 25 nVtoi W ^ JF* regions that have negative 
wavelength gravity andtopography foUow the surface topography m sign at all 

deformed terrain, on Venus. „ that of time-dependent crustal flow generates 

In contrast to the regime of tle downwelling, the topography is negative m the ear y 

topography that changes sign. 9 ve f a P^ t 1 l r Staees Q f deformation; the converse holds overa 
stagS of deformation and positive in the topography goes through zero, the geoid does 
mantle upwelling. bunng the transihon^l ^ FJJ . this cJmmhmty is not observed on Venus, 
not go tlnough zero, and the of tS geoid anomaly mimics that of the 

hi the absence of a mechanism by wh^ to ^ the regime of rapid crustal flow is 

ft criven uDwelling or downwei 1 g, y . 4 l t *!,<» inwer mantle be 


more 


nre^e.nce of large shield volcanoes in cenaui ^ regime D f rapid 

EaS are best explained as areas of mantle ^® Il S®LVerthe partial melting and ensuing 
cLstal fkm predicts crustal IhinmngovH *c“pwd outweigh the effects of crustal 

volcanism aid cmstal plutonism OTRs [l?i 

thinning on the surface topography and thus y w s Kiefer ct aL, Geophys. Ra. Lot., 

References. [11 S. C Solomon et al.. Science, 252, 2 , . 1 1 R } phaiips et aL, Science, 252, 651, 

£ f.?m^ P E^h^anet: ScLLett, &^£.’,Phy?. £* 

Phillips, J. Geophys. Res., 96, 8305, 

1992; [15] R. 1- Phillips et al., Science, 252, 651, 1991. 



am AtY) 800 0 200 4-00 600 

200 400 600 800 u Rorizontal Distance , km 

Horizontal Distance, km Th^boxis 

Figure 1. Temperature contours for the vertS^T^'and bottom temperatures are 




N94- 31110 


Cold - t. - Topograph, Ratios »" Venus-. A Globa! 
Perspective 

ttW MA 

02139) . f Terrestrial Magnetism, Carnegie 

Recently av^able^heri^ 

topography of Venus are of suf X Questions concerning the 

order 60) that they can be used to Sstan, gional scale. We 
relationship between geoid and topograpn, ^ d . to . to p 0g niphy 

have approached this question y Pj? p given a point on the 
S (6TR) on a systematic g^b J ^ For 

surface, we consider the geoi located within a specified 

gridded representation of those fie paired values we 

distance ofthe referee , PXien"eW,Lg snaightline 
determine thecoff^uon coeff ce ^ fonncr is a measure of 

The latter is the GTR at that pos * . ^^m-e j s then repeated 

the significance ofthe . de ^^^ e lding maps of the GTR and the 
for all points on the giobal^d, ^eld g the GTR on Venus, 

correlation coefficient. Un P objective and systematic search 
this approach permits us to niak J that do not 

for regions with anomalous GTOs as well as. d ^ topography . 

demonstrate any str<^ corre harmonic models ofthe 


JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 97, NO. EIO, PAGES 16.121-16,148, OCTOBER 25, 1992 


Gravitational Spreading of High Terrain in Ishtar Terra, Venus 

Suzanne E Smrekar and Sean C Solomon 

Department cf Earth, Atmospheric , and Planetary Sciences, Massachusetts Institute cf Technology, Cambridge 


Magellan altimetry measurements indicate that the mountain belts and plateau margins of Ishtar Terra 
have topographic slopes in the range 2-30°. Magellan radar images reveal that numerous sets of narrow, 
closely Spaced troughs and lineations, which we interpret as graben and normal faults, are located in 
many of these regions of high slope. The orientation of many of these graben sets perpendicular to the 
down slope direction suggests that they formed as a result of gravitational spreading. In several of the 
mountain belts, the exiensional structures are parallel to the apparent direction of shortening, consistent 
with the interpretation that gravitational spreading occurred while mountain building was still active. To 
explore the implications of this spreading for the tectonic evolution of Ishtar Terra, we model the 
process by means of a finite element algorithm for viscoelastic deformation in a vertical section of the 
crust at the margin of a plateau or broad mountain range. Flow in the crust is assumed to be governed 
by a nonlinear, depth -dependent rheology. We predict brittle failure and relaxation of topographic 
relief as functions of time, for ranges of crustal thickness (10-30 km) and thermal gradient (5-25 K/km) 
and for observed ranges in topographic elevation (1-6 km) and slope (1-30°). This parameter range 
predicts a large spread in the time scales for initial failure and significant relaxation of relief. For a 
crustal thickness greater than 10 km, a thermal gradient of 15 K/km or more, and average values of 
relief (3 km) and slope (3°), the topography relaxes to 75% of its original height in less than 10 m.y. 
Although initial failure, in the form of normal faulting on the highlands and margin and occasional 
shallow thrusting in the lowlands, is predicted to occur much earlier, values of horizontal surface strain 
that are likely to be observable (~1%) do not accumulate until significant relaxation of relief begins. 
Given lhe predicted rate of topographic relaxation, either the crust of Ishtar Terra must be anomalous in 
some way (either much stronger than predicted by standard flow laws or very thin or cold) or the 
topographic relief and high slopes have been actively built and maintained until times significantly 
younger than the average crater retention age for Venus of 500 Ma. The large apparent depth of 
compensation, the abundance of volcanism, and the uncertainty in the local age favor the view that 
mountain building in Ishtar Terra has occurred until limes at least as recent as 10 Ma. 


Introduction 

The high terrain on Venus should be subject to significant 
viscous relaxation on geological time scales [ Weertman , 
1979] because of the planet’s high surface temperature 
(~500°C) and low rate of erosion [Arvidson et al., 1991, this 
issue]. The highest topography on Venus occurs in the Ishtar 
Terra region. At the center of Ishtar Terra is Lakshmi Planum, 
a plateau with elevations 3-4 km above the mean planetary 
radius. The plateau is surrounded by the mountain belts 
Danu, Akna, Freyja, and Maxwell montes, which rise an 
additional 13, 3, 3, and 7 km, respectively, above Lakshmi 
Planum. Maintenance of these extreme elevations may be 
understood if the topographic relief of Ishtar Terra is 
geologically young or if the crust is very strong. However, 
neither of these conditions are clearly satisfied at Ishtar Terra. 
The average crater retention age for Venus, calculated from 
the density of impact craters seen in Magellan radar images, is 
approximately 500 Ma [Phillips et a/., 1991a, this issue; 
Schaber et a/., this issue]; on the basis of crater density, 
Ishtar Terra is not demonstrably younger than the average 
age. The large value of apparent depth of compensation 
(180 ± 20 km) for Ishtar Terra indicates some component of 
dynamic compensation [Grimm and Phillips , 1991], but this 
figure describes the compensation of features with 
wavelengths of at least 2000 km. Although a dynamic 
compensation mechanism at depth might be expected to 
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transmit stresses to the surface and thus contribute to recent 
changes in topographic relief, the most abrupt changes in 
elevation in Ishtar Terra occur over distances of tens to 
hundreds of kilometers and probably do not directly reflect 
scales of mantle processes. The abundance of volcanism in 
Ishtar Terra [Head, 1990; Head et a/., 1991] and particularly 
in Lakshmi Planum [Roberts and Head , 1990a] implies that 
anomalously low heat flow, and thus relatively high crustal 
and lithospheric strength, is unlikely. 

Magellan radar images (120-210 m resolution) and 
altimetry data (15-25 km footprint diameter) provide new 
insights on the stability of the areas of high relief in Ishtar 
Terra. Numerous sets of lineations interpreted as narrow 
graben, 03-10 km wide, and normal faults are visible in the 
radar images of the slopes of mountains and on the scarps 
bordering Lakshmi Planum. Much of the apparent extension 
appears to be related to topographic slope and elevation, 
indicating that it is probably the result of gravitational 
spreading. Magellan altimetry data [Pettengill et al , 1991] 
show that topographic slope is even higher than indicated 
by previous data sets [Pettengill et al ^ 1980; Barsukov et al 
1986], However, western Maxwell Montes have slopes in 
excess of 30° yet show no clear evidence of extension. 
Extreme slopes such as this are very difficult to support 
unless tectonic forces are actively building the topography. 

Even on Earth, where erosion rates are high and surface 
temperatures are much lower than on Venus, gravitational 
spreading plays an important role in the evolution of high 
topography [Dewey, 1988]. The high mountains and 
plateaus associated with the Himalayas and the Andes are 
prime examples of active gravitational spreading 
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[Froidevaux and Ricard, 1987; Dewey. 1988]. Extensions 
features are observed within and on the margins of the 
plateaus [Sebrier el al., 1985; Burchfiel and Royden, 1985]. 

In Tibet, which has an average elevation of 5 km, active 
extension occurs along a set of large (75 by 15 km), north 
trending graben in the interior of the plateau. Rom structural 
evidence, the rate of extension is estimated to be 1% m.y.* 1 
over a distance of 1 100 km, corresponding to a rate of east- 
west spreading of 10 mm/yr [Armijo el a!., 1986]. In addition 
to this Quaternary extension, geologic mapping has revealed 
fault structures indicating Miocene extension parallel to the 
direction of shortening [Burchfiel and Royden, 1985]. These 
east trending normal faults extend 800 km along strike and 
may have experienced as much as several tens of kilometers 
of displacement [Burchfiel and Royden, 1985]. North-south 
extension has been documented in the Andean Altiplano- 
Puna Plateau [Sebrier el al., 1985], which is approximately 4 
km high. Rates of Quaternary extension [Sebrier et al., 
1985] are comparable to those in the Himalayas. 

Recent studies of the formation and evolution of 
mountains have shown that gravitational spreading can be 
an important driving force for extension. Modeling of areal 
deformation is made tractable with a quasi three-dimensional, 
finite element representation of a thin sheet with vertically 
averaged rheology. England and Houseman [1989] 
employed this method to examine the transition from 
th orfning to extension as a function of thermal conditions 
«mt topographic relief in an active orogen. Bird [1989] also 
took this approach, but he included two sheets, one for the 
crust and one for the mantle, which allows for detachment 
between the layers. Vertical plane-strain models (numerical 
and analytical) have been used to model more realistic 
rheologies or specific aspects of the process of gravitational 
spreading. Bird [1991] and Kruse el al. [1991] have looked 
specifically at the role of flow of the lower crust in balancing 
stress arising from relief both of the surface and of the crust- 
mantle boundary within extensional environments. Braun 
and Beaumont [1989] combined elastic, viscoelastic, and 
elastoplastic rheologies to model the tradeoffs between the 
timing and magnitude of gravitational spreading and regional 
vertical and horizontal stresses. Most of these studies have 
tried to match the thermal and tectonic histories of specific 
rites of presently high topography, such as the Andes 
and Himalayas [England and Houseman. 1989], or where 
high topography in the past is thought to have played an 
Important role in the present tectonics, such as the Basin and 
Range Province [Braun and Beaumont, 1989; Lynch and 
Morgan, 1990; Kruse et al., 1991]. 

Modeling gravitational spreading on Venus requires a 
some what different approach. The primary constraints are 
locations of faults identified in radar images relative to the 
topography, along with broad constraints on the rheology of 
the crust and lithosphere. Regional tectonic history on 
Venus can be inferred to some degree from fault patterns 
observed in radar images [e.g.. Head, 1990; Solomon et al., 
1991, this issue; Kaula et al., this issue], but strain rates and 
thermal histories are unknown. Thus the goal of modeling is 
not to constrain the details of the thermal and tectonic 
history but rather the rheology and state of stress. 
Specifically, we address the question of how long 
topography remains after relief-building processes have 
ceased. 


In this study we model gravitational spreading of high 
topography using a plane-strain, viscoelastic finite element 
algorithm [Melosh and Rafesky, 1980, 1983]. This approach 
permits the prediction of stress, failure, and strain as a 
function of time. In Ishtar Terra, extensiona’ faulting is 
observed to occur both in the interior of the plateau and on 
the slopes of the plateau and surrounding mountains. We are 
primarily interested in deformation near the slopes, rather 
than on the interior of the plateau, for several reasons. First, 
on the margins of Ishtar Terra, extension is typically the most 
recent stage of deformation, whereas extensional features 
near the center of the plateau are often overlain by volcanic 
deposits or cut by later tectonic features (see discussion 
below). Normal faulting near the slopes is thus more likely to 
be related to present topography. Second, the deformation 
of the plateau as a whole is probably more strongly affected 
by long-wavelength dynamic compensation, which is poorly 
understood [Bindschadler and Parmentier, 1990; Grimm 
and Phillips, 1991; Lenardic et al., 1991]. Third, including 
only a portion of the plateau allows for sufficient resolution 
to predict surface strain on a scale of 1-2 km. As the 
locations of normal faults and their relationship to the 
topographic slope provide major constraints for these models, 
predicting surface strain is an important goal. 

In this paper, it is assumed that flow caused by 
gravitational spreading of high topography occurs primarily 
in the lower crust and is not strongly influenced by the 
mantle. This approach is similar to that taken by Kruse et al. 
[1991] in modeling extension in the western United States. 
Although the strain is greatest in the lower crust, we use the 
term gravitational spreading to describe both viscous flow in 
the ductile portion of the crust and the related deformation at 
the surface. Fbr the rheology of the crust, we adopt a flow 
law for diabase [Carlstan, 1982], and in contrast to the 
models cited above, we take the viscosity to be depth- 
dependent. Depth-dependent (or temperature-dependent) 
rheology is important if the effects of a detachment, or a 
ductile dlcollement, are to be included in the model. The 
stresses created by the process of gravitational spreading are 
typically small relative to most other tectonic stresses, and 
some type of weakness must be present to allow significant 
deformation at the surface. In crystalline rocks, high 
temperatures at depth in the crust can produce regions that 
are sufficiently ductile to allow gravitational spreading. 
Given the high surface temperature on Venus, gravitational 
spreading above a ductile d<collement is likely over a wider 
range of thermal gradients and crustal thicknesses than on 
Earth [Smrekar and Phillips, 1988]. While the model 
includes lower crustal flow, we do not model deformation 
involving a breakaway fault that extends from the lower 
layer to the surface and allows the upper layer to detach. On 
Earth, this type of gravity-driven deformation (sometimes 
referred to as gravitational sliding) is typically confined to 
areas where a very weak layer, such as shale, is present 
[Wiltschko, 1979]. There is no obvious mechanism for 
form in e such a weak layer on Venus. 


SETTING: ISHTAR TERRA 

Ishtar Terra (Figure 1) is unique among highland regions 
on Venus. No other highland area contains a counterpart to 
the volcanic plateau Lakshmi Planum. The mountain belts 
Danu, Akna. Freyja, and Maxwell montes surrounding 
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Fig. 1. Simplified Magellan altimetry for Iihur Ten*. The contour 
interval it 1 km. The datum for all altimetry U a radius of 6051.9 km. 
One degree of latitude is approximately 105 km. Boxes indicate the 
locations of the radar images shown in Figures 2-5 and 7-11. 

Lakshmi Planum are not found elsewhere on Venus. The 
highest point on Venus lies within Maxwell Montes, 
reaching an elevation of over 10 km above the mean 
planetary radius [Petlengill et al., 1991]. On the periphery of 
Ishtar Terra, outward from each mountain belt, lie complexly 
deformed regions classified as tessera terrains. Clotho, 
Atropos. Itzpapalotl, and Fortuna tesserae are adjacent to 
Danu, Alma , Freyja, and Maxwell montes, respectively. New 
Magellan data show details of the deformation in the 
mountain belts and adjacent tesserae [ Solomon et a/., 1991, 
this issue; Kaula et al ., this issue]. 

Pioneer Venus gravity data give an apparent 
compensation depth of 180 ± 20 km for the long-wavelength 
(greater than 2000 km) topography of Ishtar Terra [Grimm 
and Phillips, 1991]. This value is within the range found for 
many of the highland areas interpreted to be sites of mantle 
upwelling, or hot spots [Smrekar and Phillips, 1991]. 
However, Ishtar Terra differs in many ways from areas 
interpreted as hot spots, which typically display a broad 
topographic swell, rifting near the summit, and large volcanic 
constructs. Several other large highlands that have a 
generally plateaulike topographic profile, such as Ovda, 
Thetis, Tellus, and Pheobe regiones, have significandy smaller 
apparent depths of compensation, less than 100 km [Smrekar 
and Phillips, 1991]. 

The absolute age of the surface of Ishtar Terra and other 
areas of Venus is poorly constrained. Using an impactor flux 
derived from Phanerozoic terrains on Earth, Schaber et al. 
[1987] find an average crater retention age for the northern 
25% of the surface of Venus of 50-160 Ma. Applying a lunar 
model of impactor flux, in contrast, gives a much larger age of 
0.5-1 Ga [Ivanov et al., 1986; Basilevsky et al., 1987; 
Phillips et al., this issue]. Phillips et al. [this issue] conclude 
that the overall distribution of impact craters on Venus, as 
derived from the new Magellan radar data, cannot be 
distinguished statistically from a random distribution, making 
it impossible to determine relative dates for individual areas 
from crater density. However, Phillips et al. [this issue] 
argue that episodic resurfacing may nonetheless be 
consistent with the spatially random crater population, and 
they suggest that the local geological history may be used to 
infer whether the surface age of a given region differs from 


the average for the planeL Roberts and Head [1990a] argue 
that the abundance of volcanism in Ishtar Terra, and 
part icular ly in Lakshmi Planum, favors a resurfacing age that 
is less than the planetary average value. 

A wide range of hypotheses have been proposed for the 
origin of Ishtar Terra. Among these are mantle upwelling 
(hot spot) [Pronin, 1986; Grimm and Phillips, 1991], mantle 
downwelling [Bindschadler and Parmentier, 1990; 
Lenardic et al., 1991], and horizontal convergence [Head, 
1990; Roberts and Head, 1990a]. Each hypothesis must 
incorporate special circumstances to explain all of the 
characteristics of Ishtar Terra [Grinyn and Phillips, 1990, 
1991]. Mantle upwelling most easily explains the gravity 
signature but requires a lateral contrast in lithospheric 
strength to produce the surrounding mountain belts. Mantle 
downwelling can also satisfy the gravity constraint, given 
the presence of a ductile lower crust which allows for crustal 
thickening. Peripheral mountain belts may form if the 
thickening of the crust and lithosphere evolve into a second 
stage of deformation in which the lithosphere sinks, initiating 
subduction [Lenardic et al., 1991]. Horizontal convergence 
can produce the correct arrangement of mountains and 
plateau but does not satisfy the gravity data; further, the 
driving forces are not expl icit. 

Several variations on the three basic mechanisms have 
been proposed, including multiple episodes of convergence 
[Head, 1990] and convergence onto a preexisting tessera 
block [Roberts and Head, 1990a]. The latter hypothesis is, in 
part, motivated by tectonic features, including ridges and 
grooves on Lakshmi Planum, that are partially buried by 
volcanic flows [Head et al., 1990; Roberts and Head. 
1990 a,b]. The presence of a preexisting tessera block has 
also been suggested as a means of concentrating stress on 
the periphery, allowing concentric mountain belts to form 
during mantle upwelling [Grimm and Phillips, 1990]. 
Crustal underplating and lithospheric delamination have 
b ren suggested as possible mechanisms for the formation of 
Lakshmi Planum [Roberts and Head, 1990b], but these 
processes do not account for the mountain belts. 

F w-h of these models was proposed on the basis of radar 
imaging and altimetry data from Pioneer Venus and Venera 
15 and 16. Early analysis of the higher-resolution data from 
Magellan shows that the tectonics of Ishtar Terra are very 
complex and cannot easily be interpreted in terms of any of 
the above models [Kaula et al., this issue]. Further study of 
the details of the deformational history of Ishtar Terra along 
with new, higher-resolution gravity data scheduled to be 
obtained later in the Magellan mission should help to 
distinguish the relative importance of the various 
hypotheses. 

EVIDENCE FOR GRAVITATIONAL SPREADING 
* 

Criteria for Identification of Gravitational Spreading 

Gravitational spreading on Venus can be recognized by 
analogy to features observed on Earth. Although 
gravitational spreading may produce a variety of fault types 
[Dewey, 1988; England and Houseman , 1989], we 
concentrate on identifying extensional features and in 
particular on identifying graben. As will be shown below, in 
most cases gravitational spreading is predicted to lead to 
stresses sufficiently large to cause failure in extension but not 
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in compression. In mountain belts, which are produced by extensional stresses due to 

shortening, extensional strain is a more characteristic ^ J'^un pta ^ ^ ^ compressive stresses 

manifestaUon of gravitational spreading than either ^eanng gra * north-south convergence. Thus recent 

or shortening. Gravitational spreading MM generally JJ-g « ^ p rinci pally in the east-west direction 

result in a direction of greatest extensional strain parallel to extension in the Altiplano-Puna Plateau occurs in 

* "£2 s t^^Xnce the orientation of strain the northUouth direction, perpendicular to the direction 

resulting from gravitational spreading “Kiting the nature of a fault from radar data alone is 
convergent orogens. the extension direction appears to oe 



Hg. 2 a. R*iar image for a portion of w«l«m Vesu Rup». S~ jg-ljj ^ £££ 22 £ 

sinusoidal equal-area projection, with north to the top. The ^ a / 1991 ]. Numerous graben, pit chains, and collapse 
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often difficult. Characteristics suggestive of extension^ 
f ? i.itinp include linearity and sense of offset across the fault. 
Lineations that occur in pairs are a clearer indication of 
extension. A flat floor between paired, inward facing linear 
faults is diagnostic of a graben, but these characteristics may 
not always be present Some volcanic features also provide 
an indication of the orientation of the least compressive 
stress. Chains of volcanic craters, or pits, are magma 
withdrawal features fed by subsurface dikes that tend to be 
aligned perpendicular to the direction of least compressive 
stress [e!£, Gudmundsson, 1987]. The pit craters can be so 
closely spaced that they form a single, elongated depression 
£ VLSE a grabai except that the sides of the > scarps 
are scalloped rather than straight. Often lava channels 
from one end of the pit chains [see Head et at., 1991, 

Figure 22]. 

Estimation of Slopes and Strains 

In the areas identified as having experienced gravitational 
spreading, w'e estimate topographic slopes and horizontal 
surface strains. Slopes discussed below are calculated 
directly from the altimetry data points. Topographic profiles 
are given along great circles. Points along a profile we foimd 
by searching for altimetry observations within a 10-15 km 
horizontal distance and weighting the observations by*e 
distance from the point on the profile. The locus of a profile 
is only approximately a straight line on the radar imag , 
which are in sinusoidal projection. Contour m p 


accompanying radar images below are created by first 
tlK altimetric measurements to match the images, 
*e data onto a regular grid. Gridding and 
to smooth Steep slopes. It should also be 
that in regions of high slope the altimetry may mdude 
elrora due to incorrect identification of the nadir return 
among multiple returns from different points on the top. base, 

•»“ 

approximate depths of the observed graben. The lip angle 
S tiK faults making up the graben walls is assumed to be 60 
±15" On the basis of this assumption and the measured 
widths of graben walls in the radar observations. ®wt of Jte 
graben diSssed in this paper have estimated 
frv) m. Unpaired lineations with similar widths are taken to 
be normal faults with throws equal to the depths of nearby 
t^T&ror. to Urn strain estimates reflect pritnanl, 
u ncertainty in the fault dip angle. 

Scarps 

Vesta Rupes. Vesta Rupes is a scarp, 2-4 km in elevation, 
that forms the southwestern boundary of Lakshmi Planum 
between Danu and Akna montes. Most of Lakshmi Planum 
slopes gently to the south, but it rises 0.2 -1.5 km near its 
southern boundary, creating a rim along the plateau 
The character of deformation changes along the length of 
Vesta Rupes. . To the west of about 320°E, Vesta ^“pcsis 
less steep and is cut by numerous, closely spacedgraben 
co llap se features with widths of up to 8 km (Figure 2a). 
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Kg. 3a. Radar image for western Dsnu Montes and eastern VesU Rupes. The radar-bright regions near the top are part of the crest of Danu 
Montes. The plain of Lakshmi Planum appear* at the top right Near the center and lower right of the image are several long, bright, 
northwest trending lineations that are interpreted to be normal faults. Numerous, narrow, northwest trending paired bright and dark 
lineations, Interpreted lo be graben, appear in the central and upper left portion of pie Image. Both sets of lineations are approximately 
parallel to the trend of the topographic contours shown in Figure 3 b. 



330 ° 331 ° 332 ° 

fig. 3b. Contoured altimetry for the portion of western Danu Montes 
and eastern Vesta Rupes shown in Figure 3a. 


Almost all of the graben and collapse features strike 
approximately N60°W, which is similar to the regional trend 
of Vesta Rupes in this area. In the western portion of Figure 
2a, the graben parallel the scarp; in the eastern portion, the 
scarp trends more nearly east- west, and most of the graben 


intersect the scarp at an acute angle. The estimated 
horizontal strain in this area is 2-16%. Numerous pits, 
collapse features, and lava channels imply that volcanism 
accompanied extension. In many cases, the extensional 
features cut each other, often at small angles, suggesting that 
the stress field evolved over time. This region was postulated 
to be a site of gravitational sliding on the basis of Venera 15 
and 16 radar data by Bindschadler and Head [1991]; they 
also suggested that congressional ridges located 200-600 
km to the southwest of Vesta Rupes were formed as a result 
downslope compression associated with this sliding. 

To the east of 321 °E, Vesta Rupes drops more steeply to 
the plains. The easternmost portion of Vesta Rupes is shown 
in Figure 3. In this region, numerous, narrow (-1 km wide) 
graben trend approximately perpendicular to the downslope 
direction. In the extended area, topographic slope is 3.5° 
over 60 km, and horizontal strain is estimated to be 1-6%. 
These narrow faults may represent the surface expression of 
deformation above a detachment. Pit chains, isolated pits, 
and dark regions which may be flows indicate the presence 
of volcanism in the area. Some of the narrowly spaced faults 
and graben cut the radar-bright crest of the western arm of 
Danu Montes (Figure 3a), suggesting that extension 
postdates mountain building in this area. In the lower right 
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of Figure 3, several long (75 km), bright lineations appear to 
be scarps. These scarps are interpreted as normal faults, as 
they parallel graben to the west and are collinear with pit 
rfiains just to the east of Figure 3. Farther west, there are 
longer scarps (up to 400 km), parallel to Vesta Rupes, that 
appear similar to those visible in Figure 3. These features 
may also be normal faults, possibly due to gravitational 
spreading. 

Southeastern scarp of Lakshmi Planum. The edge of 
I jicshmi Planum between Danu and Maxwell Montes drops 
2-4 km to the plains and has two curved segments, convex 
toward the plateau (Figure 1). Thus the edge of the plateau 
fjin be approximately divided into northeast trending 
sections and northwest trending sections. Slopes on the 
northwest trending edges are typically steeper (8.3° over 20 
km) than those along the northeast trending edges (4.0° over 
45 km), although near 350-35 1°E there are apparent 
elevation changes of 2 km between adjacent footprints (14° 
over 8 km). Figure 4 shows a region where two such trends 
intersect, near 62°N, 350°E. North to northwest trending 
graben are perpendicular to the downslope direction in this 
region and appear to be due to gravitational spreading. 


Horizontal strain across the northwest trending graben is 
estimated to be 1-16%. Numerous pit chains and linear 
collapse features are interspersed with the graben. These 
graben occur on the raised rim of the plateau and continue 
across the slope and onto the plains at the base of the 
northwest trending scarp (Figure 4). Some of the lineations 
parallel to the northeast trending segment of the scarp at the 
eastern edge of the area shown in Figure 4 appear to be 
graben, but most of the lineations are not paired. If these 
lirv atinru: are normal faults, the estimated strain is 1-4%. The 
radar look direction is more nearly parallel to the lineations 
on the northeast trending scarp, however, which may 
contribute to their subdued appearance. Alternatively, the 
steeper slopes on the northwest trending scarps may have 
led to greater downslope extension, or a regional stress field 
may have inhibited extension on the northeast trending 
scarps while favoring extension of the northwest trending 

scarps. _ . , 

Northeastern scarp of Lakshmi Planum. The edge of 
Lakshmi Planum between Freyja and Maxwell montes 
(Figure 5) rises 0.5-1 km above the plains to the west, then 
drops 2-3 km to the east into a basin (Figure 6). A northwest 



of the plateau. 
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trending set of graben cuts a series of bright, narrowly 
spaced ridges that appear to be the result of shortening, as 
well as several volcanic flows. Both the ridges and the 
northwest trending graben are parallel to the rim of the 
plateau. The slope of the western side of the rim is 2.2° over 
15 km; slopes to the east are as steep as 13° but are typically 
4.5° over 30 km. Extension in this region is the latest stage 
of deformation and is in the downslope direction, implying a 
gravi tati onal spreading origin. On the steeper, northeastern 
slope of the plateau rim, the northwest trending graben are 
more disrupted than to the southwest, where the graben 
continue for approximately 50 km onto the plain. Horizontal 
strain in this area is estimated to be 1-5%. One large, 
northeast trending trough cuts across the eastern slope of the 
plateau edge; its trace approximately parallels the local 
topographic contour. 

Mountains 

Dam Montes. Several examples of extensional features 
are visible in Danu Montes, the least elevated of the four 
mountain ranges rimming Lakshmi Planum. One large graben 
(10 km wide, 50 km long, 0.5 km deep) occurs at a small 
nngl* to the crest of Danu Montes (Figure 7). Apparently 
compressional ridges that crosscut the southern end of this 
graben imply that the extension at the crest of the range 
occurred between episodes of compressional deformation 
and thus during the shortening that formed Danu Montes. 
To the north and west of the graben, numerous northwest 
trending collapse features and lava channels cut northeast 
trending ridges. The orientation of the collapse features 
perpendicular to the compressional structures in the 
northwestern section of the image is consistent with 
gravitational spreading within an active orogcn. 

To the north and east of the large graben, a different set of 
east-west trending pit chains cuts across the ridges at an 
oblique angle. These chains continue to the east, beyond the 
area shown in Figure 7, for a total length of about 125 km. 
Although they do not seem to be related to the topography 
in the area shown in Figure 7, to the east they are 
perpendicular to the downslope direction. The collapse 
features, lava channels, pit chains, and numerous isolated pits 
provide evidence of multiple episodes of magmatism in the 
region. 



Fig. 4c. Altimetric profile across the southeastern edge of Lakshmi 
Planum. Location is shown by arrows in Figure 4b. 


Maxwell Montes. Maxwell Montes contain the highest 
terrain and steepest slopes on Venus, with a maximum 
elevation of more than 10 km above mean planetary radius. 
Interestingly, no graben are obvious on the western slope of 
Maxwell Montes, where the local slope is as much as 30° 
over horizontal scales of several tens of kilometers. If 
unpaired normal faults are present, they are probably 
interspersed with other faults and folds, making them difficult 
to identify. Pronin and Kreslavsky [1992] interpret many of 
the topographic troughs on the west side of Maxwell 
Montes to be graben. Although this is a plausible 
interpretation, we do not include these troughs in our 
description of extensional features because none of them 
have all the characteristics of graben. 

Two perpendicular sets of graben occur on the southern 
slope of Maxwell Montes (Figure 8). The topographic slope 
in this region is approximately 2°; the elevation drops 6 km 
to the southeast over a horizontal distance of 350 km. One 
set of graben is parallel to the downslope direction; the other 
set is perpendicular to the downslope direction and curves 
along the southern margin of Maxwell Montes, giving the 
two sets a radial and concentric pattern. The graben set that 
trends parallel to the downslope direction is more prominent 
in the southern, central portion of Figure 8. Horizontal strain 
is estimated to be 1-17% across the downslope-parallel 
graben set and 1-9% across the set perpendicular to 
downslope direction. Northwest trending ridges, interpreted 
as products of shortening, are also prominent in soudiem 
Maxwell Monies. Disrupted, apparently low-relief ridges 
occur downslope from and approximately parallel to the 
northeast trending, concentric graben set These ridges, also 
interpreted as due to compression, are crosscut by the 
northwest trending graben. 

The concentric set of graben and normal faults, those 
perpendicular to the downslope direction, are probably a 
result of gravitational spreading. In many cases the graben 
cut the northeast trending ridges; in other regions both 
graben and ridges are disrupted such that it is difficult to 
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are approidlLely perpendicular to the topographic . lope and cut compression* ridges and volcanic flows. A 
bribed by Robert, and Head 11990 a) appeal, in the lower left comer. 


infer any crosscutting relationship. Patches of both radar- 
bright and radar-dark material fill areas between the ridges 
and graben. These regions of locally uniform radar 
backscatter may be volcanic flows or erosional deposits. In 
much of Maxwell Montes, as in most of the highlands (see 
also Figures 3 a, 7, 10. and 11a), the emissivity is extremely 


low [Tyler el al., 1991; Klose el al., this issue]. Thus the 
radar-bright patches could be areas that are bright because of 
a high intrinsic reflectivity. 

Two sets of narrow graben also occur on the northern 
flank of Maxwell Montes (Figure 9). A northeast trending 
set of graben trends perpendicular to the downslope 
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Rg. 6a. Contoured altimetry for southeastern Freyja Montes and 
northeastern Lakshmi Planum, including the areas shown in Rgurc 5 
(southern box) and Figure 10 (northern box). 

direction, suggesting a gravitational spreading origin; a 
northwest trending set parallels the downslope direction. 
The age relationship between the two graben sets is unclear. 
The average slope in this region is 3.8° over a horizontal 
distance of 80 kin, with a drop of 6 km in elevation. The 
horizontal strain due to extension is estimated to be 2-22% 
across the northeast trending graben set and 1-6% across the 
northwest trending set 

Freyja Montes . Numerous northwest trending graben cut 
across the southern slope of Freyja Montes, near the 
transitional region which divides die east-west trending 
ranges in central Freyja Montes from the north-south 
trending ranges in eastern Freyja Montes. The variably 
spaced graben trend generally perpendicular to the 3.9° 
regional topographic slope (measured over 40 km) and 
continue for approximately 80 km onto the more gently 
sloped volcanic deposits of Lakshmi Planum (Figure 10). 
This orientation suggests that the graben likely formed as a 
result of gravitational spreading. The estimated horizontal 
strain is 1-6%. These graben cut across the ranges of eastern 
Freyja Montes and some of the east-west ranges in southern 
central Freyja Montes. The graben are the latest stage of 
deformation in this area and cut across lineations interpreted 
to be folds and thrust faults, as well as an impact crater 
(Figure 10a). Localized volcanism is not evident in this 
region, but the smooth, dark plain of Lakshmi Planum is 



Distance (km) 

Fig. 6b. Alti metric profile across the region shown in Figure 5. 
Location is given by arrows outside the southern box in Figure 6a. 



Fig. 6c. Altimetxic profile across the region shown in Figure 10. 
Location is given by arrows outside the northern box in Figure 6a. 


generally interpreted as a series of volcanic deposits [Roberts 
and Head, 1990 ajb]. 

Extension in central Freyja Montes is more limited in 
distribution. A set of narrow, closely spaced graben trend 
north to northwest [Solomon et al. 9 1991, Figure 14]. The 
graben occur near the highest elevation in Freyja Montes, 
but they do not appear to be consistently oriented with 
respect to the local topographic slope. Typically, they are 
neither perpendicular to the downslope direction nor parallel 
to the direction of apparent shortening. In one region, where 
the graben trend northward, they are perpendicular to radar- 
light and -dark bands that have been interpreted as 
compressional ridges [ Crumpler et al 9 1986; Head % 1990]. 
Their truncation at long lineations and their curved patterns 
suggest that they have been shortened and probably 
horizontally sheared [Solomon et al. t 1991]. If these graben 
were formed as a result of gravitational spreading, they have 
been subsequently deformed such that they are no longer in 
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H« 7. Radar image of an area about 130 km wide along the northeastern arm of Danu Monte*. The highest terrain i* in the in 

the southwest comer of the image. A large graben is located in the center of the image, trendtng at a small angle to the compressional 
ridgw^Northwest trending collapse features, lava channels, and pit chains cut approximately orthogonally across the ndges. East-west 
bending pit chains also cut across the ridges. 
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low emifsivity relative to the rest of the planet (see text). 


their original orientations. No evidence of volcanism is 
present in this region. 

In eastern Freyja Montes there is abundant evidence for 
extension. Two intersecting graben sets, one trending 
northwest, the other trending north-south, cover an area 
approximately 100 km in width by 200 km in length 
[Solomon et al., 1991, Figure 16]. The northwest trending set 
appears to be the younger of the two. Local topographic 
slopes are quite high in portions of this region, up to 1 1 , 


suggesting that gravitational spreading is likely. However, 
while some of the graben are perpendicular to the downslope 
direction, others cut across the slope at a variety of angles. 

Akna Montes. Extensional faulting in Akna Montes 
occurs primarily on the western slope. Northwest trending 
graben and normal faults cut across the northeast trending 
ridges that have been interpreted as folds and faults 
[Crumpler et al., 1986; Head et al., 1990] composing the 
mountain range (Figure 1 Id). In some cases, the graben are 
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perpendicular to the folds, such as in the northern and 
southernmost portions of Figure 1 1 a; in the central region 
they cut across the trend of the folds at an acute angle. In 
the central region, the graben appear curved, due either to 
undulations in the topography or to deformation after their 
formation. If these graben are the result of gravitational 
spreading, their orientation with respect to the topography 
and the ridges indicates that they formed while the 
shortening that created Akna Montes was still active. 
Estimated horizontal strain is 1-14%. Only one small pit is 
visible in Figure 1 la, but many of the ridges in the west 
appear to be embayed, possibly by volcanic flows. South of 
the area shown in Figure 1 la, at the southern tip of Akna 
Montes, several large troughs trend approximately 
perpendicular to the downslope direction, suggesting that 
they are due to gravitational spreading. 

Lakshmi Planum 

Lakshmi Planum is dominated by volcanic flows and 
faulting, much of which is associated with the calderas 
Colette and Sacajawea. Grooves interpreted as graben are 
common [Roberts and Head , 1990a]. In some areas only 
graben occur, but in other areas, ridges interpreted to be 
compressional features are also present The age relationship 
between ridges and graben in these areas is unclear [Roberts 
and Head f 1990a]. Both ridges and graben appear to be 
closely linked to the formation of the surrounding orogenic 
belts: compressional ridges are often parallel to and graben 
are typically perpendicular to ridges in adjacent mountains 
[Roberts and Head, 1990a]. Most of the deformed regions 
are somewhat more elevated than the surrounding plains and 
thus are only partially flooded by volcanic flows [Roberts 
and Head, 1990 a]. A portion of the grooved terrain is visible 
in the southeastern comer of Figure 5. In this area, one of 
the two graben sets in the grooved terrain is parallel to the 
edge of the plateau and to the set of graben farther to the 



Fig. 8c. Altimetric profile across southeastern Maxwell Montes. 
Location shown by arrows in Figure 8 b. 


east that we have interpreted to have formed by gravitational 
spreading. Numerous, large (100 km long by 5 km wide), 
northeast trending graben intersect the southeastern scarp of 
Lakshmi Planum, hi this region, the graben are not flooded 
and appear to be the youngest deformational features, 
sometimes cutting older areas of flooded grooves and ridges. 
Gravitational spreading can easily lead to extensional 
faulting in the interior of a plateau such as Lakshmi Planum 
(see below) and may be responsible for many of the 
extensional features observed there. The orientation of the 
many of the graben in the interior of the plateau suggests 
that if the extension is the result of gravitational spreading, 
then the spreading occurred during the formation of the 
adjacent mountain belts. 

Models 

To model the time-dependent stresses and displacements 
associated with the gravitational spreading described above, 
we employ the Finite element program TECTON [Melosh and 
Rafesky , 1980, 1983]. A two-dimensional, plain-strain 
representation of the problem is assumed. An example of the 
initial finite element grid depicting a vertical section of the 
crust is shown in Figure 12a. The reference crustal thickness, 
that of the plains, is shown on the right; on the left is a region 
of elevated topography that we will refer to as the plateau, 
although it could represent the topography of either a 
plateau or a broad mountain range. A sloped region, or scarp, 
connects the plateau and the plains. The boundary 
conditions are zero vertical and horizontal velocity on the 
bottom of the grid, a free surface top boundary, and zero 
horizontal velocity and free slip vertically at the sides. The 
elements are 1 km high and 2 km wide. We adopt a 
viscoelastic rheology, a Poisson's ratio of 0.25, and a Young's 
modulus of 6xl0 10 Pa. The viscosity is calculated as a 
function of depth from a given flow law (see below) in which 
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strain rate varies approximately as a power of the deviatonc 
stress and depends exponentially on the reciprotral 
temperature. The surface temperature is taken to be 740 K, 
and a linear thermal gradient is assumed. Each row of 
elements has the same viscosity, which is eq uiv ^nt to 
assuming that the temperature is constant along both the 
bottom Mid top of the grid. The thermal gradien is thus less 
in the plateau than in the plains, by a factor equal to die ratio 
of the crustal thicknesses. The bottom boundary condmon is 
intended to simulate approximately the effects of an upper 
mantle layer that is considerably stronger than the lower 
crust As an alternative to a lower boundary that is fixed in 
both the vertical and horizontal directions, we have earned 
out calculations for a lower boundary fixed in the vertical 
direction but free horizontally; such a boundary condiuon is 


intended to simulate a detachment surface. An actual 
detachment surface is likely to display behavior somewhere 
between that of the two boundary conditions. 

The choice of a flow law for crustal material is critical to 
the estimation of time scales for deformation by gravitational 
spreading. The surface material at the Venera and Vega 
landing sites in the equatorial plains [Surkov el ai, 198 , 
19871 has a composition similar to that of basalt or dia , 
and many previous studies of crustal deformation °n Venus 
r Grimm and Solomon , 1988; Smrekar and Phillips, 1988, 
Zuber and Parmenlier, 1990; Phillips, 1990; Lenardicetal., 
1991] have adopted the flow law of Maryland diabase 
obtained from the laboratory measurements of Cjiristan 
[1982]. Some uncertainty, of course, is always introduced by 
the extrapolation of creep rates measured in the laboratory to 
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flow at geologic strain rates. Further, strain rates in rock can 
be strongly affected by the presence of water. Because of 
the extremely low water content of the lower atmosphere of 
Venus, the Venusian crust is thought to be quite dry [McGill 
etaU 1983]. Although the samples used by Caristan [1982] 
were first dried by heating, the retained water content may 

still be greater than that of the crust on Venus. The presence 
of partial melt as well as possible microcracking in the 
laboratory samples used to measure creep rates may have 
also the strength to be underestimated [Caristan, 

1982; Evans et a!., 1990]. The flow law of Caristan [1982] 
should thus probably be considered as providing only a 
lower bound on the strength of the Venusian crust As an 
approximate upper bound on the strength of the Venusian 
crust we adopt the websterite flow law of Avi Lallement 

119781 - . ^ . , 

We examine stresses and deformation over periods of up 
to several billion years. Model time steps are scaled by the 
Maxwell time, which is given by 

= cQKm / TEA 0 <j<*-» (1) 

where Q is the activation energy, R is the gas constant T is 
the absolute temperature, E is the Young's modulus, A 0 is a 
flow law constant, and O is the second invariant of the 
deviatoric stress tensor. For diabase, A 0 — 6.12x10"^ MPa^*s 
10 = 276 kJ/mole, and n = 3.05 [Caristan, 1982]. The 
parameters for websterite are A 0 - 3.16x10^ MPa*" s , Q * 
323 kJ/inole, n = 4.3 (from the data of Avi Lallement [1978] 
refit by Carter and Tsenn [1987]). The initial time step is 
tairgn to be 0.1 T«; after 10 time steps at 0.1 1«, 10 more tune 
steps are run at 1 t m , and so on, typically for 100 time steps. 
An appropriate time step for the calculations must be small 
enough to give accurate numerical results. We therefore 
adopt a value of T u given by the basal temperature and a 
stress C of 5x10® Pa. Actual values of x« in the model vary 



Pig. 9c. AMmetxic profile across northwestern Maxwell Montes. 
Location ihown by anowi in Figure 9b. 


with stress across the grid and with time, increasing by 
several orders of magnitude as the stresses relax. The 
effective viscosity also varies in timeimd space; stresses from 
the previous time step (or the initial elastic Response in the 
case of the first time step) are used to estimate the effective 
viscosity at each time* 

The finite element model was tested by comparing the 
predicted stresses and velocities to those from an analytic 
solution for one-dimensional (infinitely long, inclined slab) 
gravitational spreading [Smrekar and Phillips, 1988]. For 
finite fi^wwnt simulations with both Newtonian and power 
law, temperature-dependent rheologies, a lower boundary 
fixed in the vertical and horizontal directions, and side 
boundaries free in the vertical and horizontal directions, the 
predictions of the finite element models away from the sides 
of the grid agree very well with the analytic solutions. 

The Mohr-Coulomb criterion for basalt [Handin, 1966] is 
used to predict failure within each element The orientation 
of the principal stresses with respect to the surface 
determines the type of faulting [Anderson, 1951]. Stresses in 
our models are not adjusted in response to predicted failure. 
The horizontal strain is found by comparing element 
Hi m^ncinns before and after deformation. 

The effects of topographic slope, plateau height, thermal 
gradient, and crustal thickness on the timing of faulting and 
the relaxation of relief are evaluated by constructing models 
with a range of values for these parameters. The ranges of 
slope (1-30°) and plateau or mountain belt height ( 1-6 km) 
correspond to those observed in Ishtar Terra. Scaling the 
heat flux from Earth to Venus using the planet mass provides 
an estimate of the average thermal gradient of 10-25 KAm 
[Kaula and Phillips, 1981; Solomon and Head, 1982; 
Phillips, 1990). Most models are run for a thermal gradient 
of 15 K/km. We include some models with a thermal gradient 
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Rg. 10. Radar image of a portion of Freyja Monies. The contoured altimetry and an altimetric profile for this area are given in Figure 6. 
The northwest trending grabcn set is approximately perpendicular to the general trend of the topographic slope and represents the latest 
stage of deformation. 
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Fig. lit. Contoured altimetry of the central portion of Akna Montec 
shown in Figure 1 la. 


of 25 K/km, as well as ones with a thermal gradient of 5 K/km 
to represent the possibility of lower than average heat flow. 
The range of crustal thicknesses (10-30 km) beneath the 
plains is consistent with studies of viscous relaxation of 
impact craters [Grimm and Solomon, 1988] and models for 
the formation of regularly spaced compressional and 
extensional tectonic features [Banerdt and Golombek, 1988; 
Zuber and Parmentier, 1990]. 

We first present two examples, and then we consider some 
general results for cases with a range of crustal thickness (10- 
30 km), thermal gradient (5-25 K/km), plateau height (1-5 
km), and scarp slope (1-30°), all under the assumption that 
the crustal flow law is that of diabase. Finally, we discuss 
results for three departures from these cases: models with a 
horizontally free lower boundary condition, intended to 
simulate a detachment; models with a horizontal velocity 
boundary condition imposed on one end to mimic a regional 
tectonic stress; and models in which the crustal flow law is 

that of websterite. _ 

The history of deformation and failure associated with one 
model of gravitational spreading is shown in Figure 12. ’Die 
model incorporates a plateau height of 3 km, a bounding 
scarp of slope of 3°, a plains crustal thickness of 10 km, and a 
plains thermal gradient of 15 K/km. The finite element grid is 
shown at three model times. The first predicted failure 
(normal faulting) occurs on the plateau and on the top of the 
slope at approximately 0.1 m.y. (Figure \2b). Failure in this 
model, even at subsequent times, is confined to the plateau 
and the top of the slope and extends to a depth of only 1 km. 
Although the greatest compressive principal stress in the 
plains is horizontal, no thrust faulting is predicted. 
Significant relaxation of the topography of the plateau, 
resulting in a decrease in the angle of the slope of the 


bounding scarp, docs not occur until 100 m.y. (Figure 12c). 
The finite element grid is extended under the plateau and 
compressed beneath the plains, a result of flow concentrated 
at the base of the crustal layer. At 100 m.y. the failure 
criterion is no longer satisfied. The horizontal surface strain 
increases with time as long as the topography is still relaxing; 
at 100 m.y. the maximum accumulated horizontal surface 
strain is approximately 9% and occurs in the plateau. 

Failure and relaxation of relief occur much more rapidly for 
models with thicker crust and a higher plateau. Figure 13 
shows the results from such a model with a 30-km-thick 
crust, a plateau or mountain belt height of 6 km, a surface 
slope of 30°, and a plains thermal gradient of 15 K/km. A 
height of 6 km (above Lakshmi Planum) and a slope of 30° 
are typical of western Maxwell Montes. The first failure on 
the plateau occurs at approximately lfr 2 yr, the slope first 
fails at 10 yr. The pattern of deformation in Figure 13 is quite 
different from that shown in Figure 12. Failure near the 
surface of the plateau does not extend to the plateau edge, 
and failure near the surface of the scarp is at the base of the 
slope rather than near the top (Figure l3 b). Thrust faulting in 
the plains occurs at distances of greater than approximately 
40 km from the base of the scarp. The plateau has relaxed to 
less than 75% of its original elevation by 0.01 m.y. (Figure 
13c). At this time normal faulting is predicted at depths of up 
to 3 1cm beneath the plateau and at depths of 1-2 km near the 
base of the scarp. Ductile flow, which occurs much more 
rapidly for thick crust than for thin crust for a given thermal 
gradient, is greatest at the base of the crust. The region of 
highest strain is localized in the bottom row of elements, with 
the greatest deformation occurring beneath the scarp. 

The timing of predicted initial normal faulting beneath the 
plateau is largely controlled by the effective viscosity at the 
bfl y of the crust. Figure 14 shows the time for the first 
failure on the plateau as a function of plateau height and 
crustal thickness. Increasing the crustal thickness from 10 to 
20 km or from 20 to 30 km causes failure to occur at times 
earlier by a factor of lO^-lO 3 . There is a large difference in 
initial failure time for plateau heights of 1 and 3 km (factor of 
10M0 5 ), but no more than a factor of 10 difference in failure 
time for heights between 3 and 5 km. Temperature at the 
base of the crust increases with crustal thickness, and 
deviatoric stress near the base of the crust increases with 
both crustal thickness and plateau elevation. The nonlinear 
change in predicted failure time is due both to the 
exponential dependence on temperature and to the power 
law dependence on deviatoric stress of the effective 
viscosity. 

Crustal thickness and plateau height have similar effects 
on the timing of initial failure on the scarp. Predicted first 
failure on the scarp always occurs at the same time as, or 
somewhat later than, failure on the plateau. Typically, for 10- 
km-thick crust, failure on the plateau and scarp are 
simultaneous. For 20- or 30-km-thick crust, the time of first 
failure on the scarp can lag behind that on the plateau by up 
to a factor of 10 3 . In a few cases, no failure occurs on the 
scarp, even though failure is predicted on the plateau. 

The effect of the slope of the scarp was also studied by 
comparing cases with a plateau height of 3 km, crustal 
thicknesses of 10, 20, and 30 km, and slopes of 1°, 3°, 10°, 
and 30°. For crust of 10 to 20 km thickness, the effect is 
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small. For a 1° slope, failure on the plateau is predicted to 
occur up to a factor of 10 later than for a slope of 3®. 
Increasing the slope up to 30° does not have a significant 
effect on initial failure times. For a 30-km-thick crust, the 
effect of scarp slope, however, is much more pronounced. 
Increasing slope from 1° to 3° results in a decrease in 
predicted initial failure time by a factor of 10 3 . Between 3° 
and 10° slope and between 10° and 30° slope, the decreases 
are factors of 10. 

Failure by thrust faulting in the plains is predicted in some 
cases with plains crustal thicknesses of 20-30 km and plateau 
heights of more than 1 km. For a given thermal gradient, 
predicted thrust faulting occurs earlier for larger crustal 
thickness and plateau height Initial thrust faulting occurs 
laer than first normal faulting in the plateau by a factor of 
10M0 4 , typically not occurring until the topography begins 
to relax significantly. In all cases, thrusting occurs at 


distances of at least 20-50 km from the base of the scarp and 
only within the upper 1 km of the plains crust The width of 
the zone of predicted thrusting is likely to be controlled by 
the width of the plains in the model, as thrusting typically 

continues to the edge of the model. 

The accumulated horizontal strain at the surface at the time 
of predicted first failure is quite low, typically 0.01% The 
horizontal strain remains small until the topographic relief 
begins to relax significantly; strain is approximately 1-25% 
when the topography has relaxed by 25%. Strain is most 
strongly affected by ciustal thickness, with thinner crust 
resulting in higher strain. Strain increases typically by a 
factor of 3-10 for a decrease in crustal thickness from 30 to 
0 km. As relaxation of the topographic relief occurs 
primarily through flow in the lower crust, with a thicker crust 
near-surface material is less coupled to the deformation at 
depth and experiences less horizontal strain. The calculated 
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strains may be underestimates of the actual strain that would 
be produced by gravitational spreading, inasmuch as the 
model does not explicitly include strain by faulting. Once a 
fault forms, deformation may proceed more readily by 
movement along the fault than by creating a new fault. 


leading to higher strain than predicted by these continuum 
models. 

The first predicted extensional failure in the model may not 
produce strains sufficiently large to be comparable to 
observed strains. Minimum estimates of horizontal strains 
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Fig. 13. (a) Initial finite element grid for a model with a plateau elevation of 6 km, a bounding scarp of dope 30°, a plains crustal thickness 
of 30 km, and a thermal gradient in the plains of 15 K/km. The vertical exaggeration is 2:1; the elements are 1 km high and 2 km wide, (b) 
Gnd ax approximately 100 yr. Thrust faulting is predicted in the plains farther than 40 km from the scan) (c) Grid at aDoroximalelv 
10,000 yr. The plateau elevation has decreased by approximately 25%. F y 
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Hg. 14. The time at which failure (by norma] faulting) first occurs on 
the plateau is shown as a function of plateau elevation for three values 
or plains cnistal thickness. The initial slope of the scarp in each 
modd is 3 , and the plains thermal gradient is 15 K/lcm. Open 
symbols represent models with a diabase flow law [Caristan, 1982 ]' 
*ohd symbols are for a websterite flow law [Avi Lallement, 1978 ], All 
times fa Figures 14-15 are given only to the nearest order of 
magnitude in seconds. 


documented in regions of extension in Ishtar Terra are 
approximately 5-10%. It is unlikely that a surface strain of 
0.01%, as is typical at the time of predicted initial failure, 
would be recognizable in the radar images as extensional 
faulting. A resolvable value of horizontal surface strain may 
be closer to 1 %. The models predict that horizontal strain at 
the surface on the order of 1% is not produced until much 
later than the time of initial failure. In most cases, 1% 
accumulated strain at the surface occurs later than initial 
failure by a factor of 10M0 5 . For a thermal gradient of 15 
K/km, 1% cumulative strain is usually attained in 10M0 8 yr 
although in some cases with a 30-km-thick crust horizontal 
surface strain never reaches 1%. As discussed above, the 
higher the effective viscosity of the crust, the higher the 
strain. Thus for many cases with a thermal gradient of 25 
K/km, 1% horizontal strain at the surface is not reached, even 
though the topography relaxes significantly. 

The timing of relaxation of topographic relief is primarily 
affected by the crustal thickness and the associated 
temperature at the base of the crust. The time at which the 
average elevation in the plateau has relaxed by at least 25% 
is shown as a function of plateau height and crustal 
thickness, for thermal gradients of 15 and 25 K/lcm, in Figure 
15. For a given thermal gradient and crustal thickness, 
predicted relaxation times for plateau heights of 1, 3, and 5 
km vary by a factor of 10 at most, with a decrease in 
relaxation time typically occurring between heights of 1 and 
3 km. Typically, a model with a crustal thickness of 20 km 
will relax faster (by a factor of 100) than a model with a 
crustal thickness of 10 km. The difference in relaxation time 
between models with crustal thicknesses of 20 and 30 km is 
also a factor of 10-100; cases with higher thermal gradient 
tend to have larger differences. Larger differences in 


relaxation time are predicted for models with the same 
plateau height and crustal thickness but with thermal 
gradients that differ by 10 K/km, as is consistent with the 
larger changes in basal temperature. For a thermal gradient 
of 25 K/km and a plateau height of 1-5 km, all the models are 
predicted to relax by at least 25% within 100 m.y. for a 10- 
km-thick crust, 0.1 m.y. for a 20-km-thick crust, and 1000 yr 
for a 30-km-thick crust A thermal gradient of 15 K/km gives 
relaxation times of up to 1 b.y., 100 m.y, and 0.1 m.y. for 
crustal thicknesses of 10. 20, and 30 km, respectively. This 
strong dependence on temperature is because the flow is 
largely controlled by the lowest viscosity in the crustal 
column, which in turn is controlled by the basal temperature. 

An interesting consequence of the relaxation of the 
topographic relief is the development of topographic profiles 
that resemble those associated with lithospheric flexure. As 
the relief in the plateau decreases, the plains at the base of 
the scarp develop a foredeep and outer rise and are predicted 
to fail in extension (Figure 13c). 

As discussed above, detachment surfaces in regions of 
gravity-driven deformation on Earth develop most often in 
layers of soft, sedimentary rocks. Detachment surfaces in 
crystalline rocks, however, are also known in such structures 
as thrust belts and metamorphic core complexes [Armstrong 
and Dick, 1974; Bhattacharya, 1987; Miller et al., 1983). 
Although not well understood, mechanisms such as strain 
softening [Schmid, 1983] are thought to be important in 
reducing the strength within the detachment zone and 
localizing deformation. To simulate a weak detachment 
surface, we use a lower boundary condition that is free in the 
horizontal direction and fixed in the vertical direction. 

Predicted failure times are strongly influenced by the 
horizontally free lower boundary condition, with failure on 
the plateau occurring at times earlier by a factor of 103-106 
than those shown in Figure 14. The difference in predicted 
initial failure time between the horizontally free and fixed 
lower boundary conditions is least for those cases with a 
thick crust and high plateau. For these cases, the low 
viscosity region at the base allows deformation to occur 
rapidly even with a fixed lower boundary. Differences in the 
time of significant relaxation of relief between models with 
different lower boundary conditions are less than for initial 
failure time. With a horizontally free lower boundary, 25% 
relaxation of the initial plateau elevation occurs either at the 
same time as or earlier by as much as a factor of 100 than for 
models with a fixed lower boundary. Thus if a weak 
detachment does develop at depth, the times for initial failure 
and relaxation shown in Figures 14 and 15 will be 
overestimates. 

Regional tectonic forces can strongly influence the time 
scale of gravitational spreading. To simulate a regional force 
opposing gravitational spreading, the boundary condition at 
the plains end of the finite element grid was modified to 
include a vertically uniform horizontal velocity directed 
toward the plateau. To avoid failure of the finite element 
code due to excessive changes in the shape of individual 
elements, we use a lower boundary condition that is free in 
the horizontal condition. Two cases were considered. The 
first included a 20-km-thick crust, a plains thermal gradient of 
15 K/km, a plateau height of 3 km, and a scarp slope of 3°; 
the second included a 30-km-thick crust, a thermal gradient 
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of 15 K/km, a plateau height of 6 km, and a scarp slope of 
30°. Initial extensional failure develops as it would if the 
velocity boundary condition were zero. In time, however, 
the shortening induced by a nonzero boundary velocity 
overcomes the extension driven by the topographic relief, 
and failure is predicted to occur only by thrust faulting. For 
both cases, a boundary velocity of 100 mm/yr will prevent 
gravitational spreading after a period of approximately 10 4 
yr. If the velocity is significantly less than this value, there 
will be little change in the gravitational spreading on 
geologic time scales from models without an imposed 
velocity boundary condition. 

In a final set of models we incorporated the websterite 
flow law of Avi Lallement [1978] and Carter and Tsenn 
[1987] to simulate a crust much stiffer than would be given 
by the flow law for diabase. For a fixed lower boundary 
condition and the same ranges of crustal thickness, thermal 
gradient, and plateau height, initial failure is predicted to 
occur later, typically by a factor of lO^-lO 5 , than in similar 
models with a diabase flow law (Figure 14); failure times for 
thin crust or a low plateau differ by much less (Figure 14). 
Times for significant relaxation of relief are also much 
different for a websterite rheology than a diabase rheology. 
For the same parameter ranges as the cases shown in Figure 
15, none of the models with the websterite flow law 
approached 25% relaxation within 1 b.y. Even with a 
thermal gradient of 25 K/km, only those models with a crustal 
thickness of 30 km reached 25% relaxation in 100 m.y. to 1 
b.y. 

Discussion 

The models presented in this study indicate that 
gravitational spreading is likely to occur in most of the high 
terrain of Ishtar Terra, on time scales less than the average 
crater retention age of Venus, for observed values of 
topographic slopes and elevations and likely values of 
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Fig. 15. The lime at which ihe height of the plateau has relaxed by 
25% is shown as a function of plateau height for three values of plains 
crustal thickness. The initial slope of the bounding scarp is 3°. The 
open symbols are for a plains thermal gradient of 15 K/km, solid 
symbols are for a thermal gradient of 25 K/km. 


crustal thickness and thermal gradient. This result implies 
that regions of high elevation either have anomalously cold, 
strong crust or were recently (or are currently) subject to 
regional tectonic forces that are greater than those driving 
gravitational spreading. In the former case, this study 
provides a lower bound on the strength of the crust. In the 
latter case, the models predict the time period over which 
high topography has been relaxing. To constrain the 
rheology or the duration of gravitational spreading, we must 
some assumptions about the age of Ishtar Terra and the 
possible degree of relaxation of the topography. As 
discussed above, the absolute age of the surface of Ishtar 
Terra is quite uncertain. As a basis for discussion, we 
consider ages of 50-500 Ma; the former is the surface age in 
the northern 25% of the planet derived from a terrestrial 
impactor flux rate [Schaber et al., 1987], and the latter is the 
average crater retention age for the planet if the rate of crater 
production is scaled from the lunar impactor flux rate 
[Phillips el al., this issue]. As the elevation of Ishtar Terra is 
much higher than for most of Venus and scarp slopes are 
generally large, we also make the assumption that the 
topography has not relaxed significantly. Specifically, we 
examine model predictions for 25% relaxation of the initial 
relief. Given that our models may underestimate surface 
strain, it is reasonable to assume that observable strain will 
accumulate by the time 25% relaxation of the topography 
occurs. 


Age Constraints 

If the crust in Ishtar Terra is not anomalously strong 
relative to predictions from laboratory measurements on 
diabase, then regional tectonic forces are likely to have acted 
to maintain or enhance topographic slope within the last 10 7 - 
10 s years. For most of the regions' where gravitational 
spreading is believed to have occurred, the relief is 
approximately 2-3 km; at Maxwell Montes, 5-6 km is a better 
approximation. For a thermal gradient of 15 to 25 K/km, only 
those cases with a crustal thickness of 10 km require times of 
the order of 10-100 m.y. to relax by 25% (Figure 15). For a 
thermal gradient of 15 K/km, those cases with a 20-km-thick 
crust relax in less than 10 m.y.; those with a 30-km-thick 
crust relax in less than 1 m.y. For a thermal gradient of 25 
K/km, cases with a 20- or 30-km-thick crust relax within 1 
m.y. Almost all cases require 100 m.y. or more to relax by 
25% for a thermal gradient of 5 K/km. Thus for thermal 
gradients of at least 15 K/km and crustal thicknesses greater 
than 10 km, our models predict 25% relaxation within several 
million years. 

Given the high topography (up to 10 km of relief at 
Maxwell Montes) and steep slopes, it is likely that some 
crustal thickening has occurred and that the crustal thickness 
exceeds 10 km in many areas of Ishtar Terra. The overlap of 
volcanism and extension in many areas of Ishtar Terra also 
suggests that the thermal gradient is unlikely to be as low as 
5 K/km. Thus the time scales for relaxation imply that 
orogeny in many areas of Ishtar Terra is likely to have 
continued until at least several million years ago. 

The presence of extensional faulting does not necessarily 
imply that regional tectonic forces are no longer active. On 
Earth, episodes of shortening and extension are believed to 
overlap in time. Examples of this process of 




SMREKAR AND SOLOMON: GRAVITATIONAL SPREADING OF ISHTAR TERRA 


16,143 


contemporaneous shortening and extension have been 
documented in the Andes and in the Himalayas. Regional 
compressional forces act to increase the topographic height 
and slope to the point where gravitational spreading forces 
are able to balance the tectonic forces. Once spreading has 
occurred, the regional compressive forces may again be 
dominant. In Ishtar Terra examples of apparently 
contemporaneous extension and shortening are found in 
Damt, Freyja, and Akna monies (see Figures 7, 10, and 11a, 
respectively). 

Rheological Constraints 

Conversely, if we assume that the age of the topography 
in I s ht a r Terra is at least 50 Ma and that no regional tectonic 
forces have opposed gravitational spreading during that time, 
then the time scales for relaxation of the topography provide 
a constraint on the Crustal thickness and strength. 
Topography in Ishtar Terra will not relax by more than 25% 
over several tens of millions of years if the temperature at the 
base of the crust is less than approximately 1000 K for a 
diabase flow law. This temperature corresponds to an 
effective viscosity at the base of the crust of approximately 
10 11 Pa rK For thermal gradients of 15 and 25 K/km, the 
crust must be 10 km or less in thickness (see Figure 15). 
Crustal thicknesses of 20-30 km will not relax in several tens 
of millions of years if the thermal gradient is as low as 5 K/km. 
As there is evidence for probable gravitational spreading in 
Ishtar Terra, a maximum time of several billion years for 
significant relaxation of relief provides an upper bound on 
the strength of the crust. For 25% relaxation to occur within 
several billion years, the temperature at the base of the crust 
must be greater than approximately 880 K, corresponding to 
an effective viscosity of 10 12 Pa s* 1 . For example, a crustal 
thickness of 10 km and a thermal gradient of 5 K/km will not 
lead to relaxation on this time scale. If the strength of the 
crust is closer to that predicted by a websterite flow law, the 
temperature at the base of the crust must be near 1400 K for 
significant relaxation to develop within several billion years. 
Of those models with a websterite flow law, only those with 
a thermal gradient of 25 K/km and a crustal thickness of 30 
km predict significant relaxation within this time scale. 

Another potential rheological constraint involves the 
accumulated horizontal surface strain predicted to develop as 
the topography relaxes. As shown above, the strain is 
sensitive to die crustal thickness and thermal gradient, which 
together control how well coupled flow in the lower crust is 
to deformation at the surface. Models with a crustal 
thickness of 30 km develop surface strains that are typically 
about 0.1-1% at the time when relaxation reaches 25%. Thus 
if the topography in a region appears to have relaxed but 
does not exhibit evidence for extensional faulting, it may 
have a very thick crust. This situation may occur in 
mountain belts with large crustal roots. However, two 
factors make application of this constraint very uncertain. 
The first is the interpretation of predictions of horizontal 
surface strain by models that do not include the effects of 
faulting. The second is the uncertainty in the estimate of 
observed surface strain. This factor does not present as 
severe a limitation, however, as strain inferred from 
observations of faulting is more likely to be an underestimate 
than an overestimate. 


Comparison With Observations 

The models of gravitational spreading discussed above 
predict the distribution of faulting, the type and depth of 
faulting, the horizontal strain, and the shape of the 
topographic profile. All of these characteristics can be 
observed or estimated from Magellan altimetry and radar 
images. The general agreement between the strains estimated 
from the images and those predicted by the models was 
discussed above. The models predict ihat regions of high 
slope (-10°) will fail first through normal faulting at the base 
of the slope, while scarps of lower slope will fail first at the 
top (Figures 12 and 13). In most of the areas of extensional 
faulting in Ishtar Terra, such simple relationships are not seen. 
However, in a few areas where extension constitutes the 
only apparent tectonic activity, the relationship predicted by 
the models appears to hold. The best example is the 
southeastern scarp of Lakshmi Planum (Figure 4) along the 
northwest trending limb of the scarp, which has an average 
slope of 8°. Downslope extensional strain is visible at the 
base of the scarp, on the scarp, and at the top of the plateau. 
Along the northeast trending limb of the scarp, the average 
slope is approximately 3°, and lineations interpreted as 
possible normal faults occur only at the top of the scarp. 

At the base of the scarps along the southern and 
southeastern margins of Lakshmi Planum there are frequent 
sets of long, narrow, anastomosing ridges (e.g„ Figures 2a 
and 4a), which we interpret to indicate crustal shortening of 
modest magnitude. This type of low-strain, compressional 
deformation at shallow depths is predicted by the models to 
occur at distances at least 20-30 km from the base of the 
scarp. Although these ridges sometimes occur somewhat 
closer to the base of the bounding scarp of Lakshmi Planum 
than this range (Figure 4a), their positions and strains are 
broadly consistent with model predictions. 

The range of predicted depths of extensional faulting is 1- 
10 km beneath the plateau, 1-3 km beneath the scarp, and 1-3 
beneath the base of the scarp. Greater plateau elevations 
and crustal thicknesses result in greater depths of faulting 
beneath the plateau. The average spacing between normal 
faults in Ishtar Terra is 1 to 15 km. A greater depth extent of 
surficial normal faulting would be expected to be reflected in 
a greater spacing between faults, although this relationship 
can be complicated by such effects as listric faulting. A 
general rule of thumb is that the depth extent of faulting is 
equal to the characteristic fault spacing divided by 2 n. This 
rule gives a depth extent of approximately 0.1 to 2 km. 
Given the uncertainties in this general rule, this range is in 
general agreement with the model predictions. 

The importance of gravitational spreading on Venus has 
previously been discounted on the basis of an expected 
failure pattern which includes a breakaway fault at the top of 
a slope, extensional faulting on the slope, and shortening at 
the base of the slope [Bindschadler and Head , 1991]. While 
such relationships are observed for gravitational spreading in 
soft sedimentary rocks on Earth [Pederson, 1987], the models 
in this study do not display this pattern of deformation. 
Ductile flow at depth does bear some resemblance to a 
ddcollement, but strain at the surface remains small. For die 
cases considered here, lower crustal flow does not lead to a 
breakaway fault, as surface strains are low and the zone of 
faulting does not extend to depths of significant flow. 
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However, in cases where flow extends over much greater 
distances, it is possible that strains would be sufficient for a 
breakaway fault to be produced. A breakaway fault has not 
yet been identified in the Magellan radar images; however, 
such a fault might be difficult to discern. 

As noted above, as the topography relaxes the plateau 
margin develops a profile resembling that of a flexural bulge. 
A good example of such profile is observed along the 
southeastern scarp of Lakshmi Planum (Figure 4c). This 
profile may reflect upwarping of the plains in response to 
partial viscous relaxation of the plateau edge. The 
extensional faulting observed in this area may be due to the 
upwarping. Profiles which resemble flexural profiles are 
found in many areas of Venus. In Ishtar Terra, examples are 
found outward of Uorsar Rupes to the north of Freyja 
Montes, at western Maxwell Montes, and at eastern Akna 
Montes [Solomon et a/., this issue]. Similar topographic 
profiles are often found exterior to coronae [Stofan et al. 9 
this issue]. Solutions for the elastic or mechanical 
lithospheric thickness in these regions sometimes greatly 
exceed expected values for Venus [Sandwcll and Schubert , 
this issue]. 

As a basis for comparison with estimates of the thickness 
of the elastic lithosphere derived elsewhere on Venus from 
topographic profiles assumed to be the result of lithospheric 
flexure, we fit the topographic profile for the plains portion 
of our models to the analytical solution for flexure resulting 
from the superposition of a bending moment applied to the 
end of a broken elastic plate and a line load on a continuous 
elastic plate [Sandwcll and Schubert , this issue]. An elastic 
plate thickness of approximately 10 km fits the profile shown 
in Figure 13c with a root-mean-squared error of 25 m. This 
value of elastic plate thickness should not be compared 
directly with values obtained by fitting actual topographic 
profiles, because our models do not include the effect of 
deformation in the mantle. Ductile flow of the lithosphere on 
Venus may be an important process for modifying 
topography, however, and can lead to profiles similar to 
those arising from elastic flexure. 

Comparison With Gravitational Spreading on Earth 

Active gravitational spreading in southern Tibet occurs in 
the form of approximately 15-km-wide, north-south trending 
graben with a typical spacing between sets of faults of 200 
km. The estimated surface strain in this 1 100-km-wide area is 
1% [Armijo et a/., 1986]. In the Andes, extension occurs 
across a 400-km-wide area, where fault spacing is 
approximately 1-20 km; surface strain is also estimated to be 
1% [Sebrier et c/., 1985]. In contrast, extension in Ishtar 
Terra is confined to areas not more than 200 km across, and 
estimated strains are greater, in the range of 5-15%. Average 
graben widths (1-10 km) and spacings (5-20 km) are 
comparable to those in the Andes but somewhat less than 
those in the Himalayas. The greater surface strains on Venus 
may be a result of a lesser crustal thickness. 

Additional Effects 

Several effects not included in this study may also 
influence the process of gravitational spreading. The surface 
temperature in the mountains is lower than the temperature 
of 740 K used in this study, as the atmospheric temperature 


decreases by 7.7 K/km with increasing altitude [Se{ff t 1983]. 
Areas considered in this study have elevations up to 8 km 
above the mean planetary radius, corresponding to a 
decrease in temperature of up to 60 K. Although we did not 
specifically incorporate this difference, the effect is mimicked 
to some degree because the temperatures at the base and at 
the surface of the crust are the same under the plateau as 
under the plains in our models, giving a lower thermal 
gradient in the plateau. If modeled explicitly, the lower 
surface temperature would probably result in slightly longer 
times for predicted initial failure and lower surface strains. 
However, as discussed above, the temperature at the base of 
the crust is the primary variable controlling the rate of 
deformation. 

The nature of the compensation of the topography is likely 
to have a major effect on gravitational spreading. In the 
models of this study, the topography is not compensated, and 
no flow occurs deeper than the base of the crust. The long- 
wavelength topography (>2000 km) of Ishtar Terra has an 
apparent depth of compensation of 180 km [Grimm and 
Phillips , 1991], but no information is available on shorter 
wavelengths. For comparison, the models in this study have 
horizontal dimensions ranging from 140 to 300 km. Such 
short-wavelength topography in Ishtar Terra may b.e 
isostatically compensated or maintained partly by regional 
tectonic stresses arising from mantle convection. Flexure is 
likely to be important at even shorter wavelengths than 
those considered here. 

Compensation of topography by Any isostasy would lead 
to shorter times of initial failure and significant relaxation of 
relief for most of the cases considered. For a crustal density 
of 3000 kg/m 3 and a mantle density of 3300 kg/m 3 , an 
elevated region is balanced by a crustal root of thickness 10 
times that of the excess elevation. Thus, for cases with a 1- 
km-high plateau, crustal compensation might occur at depths 
within the range considered (10-30 km), depending on 
crustal thickness. Higher elevations that are completely 
compensated by crustal thickness variations would have 
roots that extend to depths greater than 30 km. The higher 
temperature at the base of a deep crustal root would lead to 
more rapid relaxation, if the region is in thermal equilibrium 
with surrounding areas. 

The effect of dynamic compensation is less clear, as the 
stress field and crustal thickness evolve with time. Mantle 
upwelling or downwelling may be responsible for 
compensating the long-wavelength topography of Ishtar 
Terra and may indirectly compensate shorter-wavelength 
variations in the topography by inducing variations in stress 
that can result in differences in crustal thickness. Mantle 
downwelling beneath a weak crust can produce crustal 
thickening [Bindschadler and Par mender , 1990; Lenar die 
et al„ 1991]. Mantle upwelling is typically expected to thin 
the crust, but if lateral heterogeneities in lithospheric strength 
are present, crustal thickening over limited regions is also 
possible [Grimm and Phillips , 1991]. Generation of new 
crust by partial melting above a zone of upwelling may also 
increase the crustal thickness [Phillips et a/., 1991b]. Thus 
both upwelling and downwelling models call on crustal 
thickening to produce the observed topography in Ishtar 
Terra. Further, the inclusion of stresses due to dynamic 
processes would require a larger-scale model than considered 
in this study. 
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Kiefer and Hager [1991] present a one-dimensional model 
of gravitational spreading that includes a bottom velocity 
condition meant to simulate traction on the base of the crust 
due to mantle downwelling beneath Ishtar Terra. By 
balancing the flux out of the region due to gravitational 
spreading, for a fixed topographic slope and elevation, with 
the flux into the region specified by the bottom boundary 
condition, they find the bottom velocity that will lead to no 
nd flux. Their results show that for parameters appropriate 
for the Fortuna Tessera region and basal velocities in the 
range of 0-500 mm/yr, only anomalously thick or 
anomalously cold crust will fail to undergo gravitational 
spreading. Although their model is a steady state calculation, 
and thus not easily comparable to results from this study, it 
does suggest that gravitational spreading is possible even in 
the presence of mantle downwelling. 

Incorporating a mantle layer into the study would 
probably change the results somewhat, even if convection 
does not contribute significantly to the stress field. Kruse et 
al. [1991] considered the problem of stresses due to high 
topography in conjunction with regional extension. They 
found that allowing flow in a mantle layer to accommodate 
deformation of a crust-mantle boundary caused the crust to 
thicken beneath extensionally thinned areas at rates 30-50% 
slower (for crustal thicknesses of 10-30 km) than for the case 
where a rigid lower crustal boundary was assumed. The 
results of Kruse et al. [1991] may not be directly applicable 
to Ishtar Terra because they were obtained for constant- 
viscosity layers in an extensional stress field and for lower 
temperatures than those appropriate for the crust of Venus, 
but they do suggest that including the effect of mantle flow 
would result in somewhat lower relaxation rates than found 
in this study. 

Implications for the Tectonic History of Ishtar Terra 

This study bdicates that two quite distinct possibilities 
exist for Ishtar Terra: either the crust is anomalously strong 
and no regional tectonic forces are required to support the 
topographic relief or tectonic forces have continued to build 
the topography and steepen slopes until geologically recent 
time. Crater retention ages, gravity data, and the relationship 
between volcanism and tectonism provide independent 
information. 

As discussed above, an absolute surface age for Ishtar 
Terra is quite uncertain. Crater retention ages are meaningful 
only for areas much larger than even the large-scale tectonic 
features of Ishtar Terra [Plaut and A rv ids on, 1988; Phillips 
et al., 1991 a, this issue]. Geologic information about the 
surface age is unclear. The 105-km-diameter crater, 
Cleopatra, in Maxwell Montes is a particularly obvious 
example of an undeformed impact structure in a region of 
very high topography [Phillips et al., 1991a]. However, 
other impact structures in the region have experienced 
tectonic modification (e.g., extension across the crater visible 
m Figure 10). As surface strain associated with gravitational 
spreading is small, this type of deformation may not be very 
efficient at reducing the crater population or even 
substantially modifying crater shapes. Thus the estimates of 
the age for Ishtar Terra and the surrounding region may not 
provide a strong constraint on the recency or duration of 
gravitational spreading. 


As noted earlier, gravity data for Ishtar Terra indicate that 
the long-wavelength topography is dynamically supported 
in the mantle [ Grimm and Phillips, 1991]. Dynamic 
compensation, however, does not specifically rule out the 
possibility of cold crust. Downwelling might have led to 
thickened, initially cold crust; eventually, however, the crust 
will equilibrate to more typical temperatures. 

A low temperature at the base of the crust in Ishtar Terra 
may not be consistent with the extent and timing volcanism 
found in the region. In areas where pits, pit chains, and 
collapse features overlap with graben and extensional 
faulting, the volcanism appears to be contemporaneous with 
gravitational spreading. Of the regions of extensional 
faulting described in this study, the evidence for volcanism is 
least strong for Freyja Montes and northern and southern 
Maxwell Montes. One interpretation is that the crust is 
anomalously cool and strong because of recent shortening in 
• these regions. 

Last, we must also consider the possibility that because the 
crust on Venus may be extremely dry, the lower crust is much 
stiffer than indicated by the laboratory flow law obtained for 
diabase [Caristan, 1982]. As demonstrated by models with a 
websterite flow law, if the rheology is much stiffer, high 
elevations with moderate crustal thicknesses and thermal 
gradients might support gravitational stresses for periods of 
up to several billion years. However, this possibility raises 
the question of why Ishtar Terra is the only region of 
extremely high topography preserved over the geologic 
history of Venus. Under the assumption that gravitational 
spreading and viscous relaxation of high terrain have been 
important processes in the recent geologic history of Venus, 
Ishtar Terra represents a variation on landforms observed 
elsewhere on the planet. Gravitational spreading is also 
believed to part of the evolutionary sequence of tesserae and 
coronae [Bindschadler et al., this issue; Squyres et al., this 
issue]. One implication of a very stiff crust on Venus would 
be that the crust could be thicker than previous estimates of 
10-30 km. 

We favor the hypothesis that many of the areas of greatest 
relief and steepest slopes in Ishtar Terra have recently been 
tectonically active. The primary evidence against this 
interpretation is the estimated crater retention age for Venus. 
We note, however, that the recency of tectonic activity does 
not provide a means to distinguish among the major 
competing hypotheses for the origin of Ishtar Terra. High- 
resolution gravity data for Ishtar Terra, which are scheduled 
to be obtained later in the Magellan mission, should provide 
further constraints on plausible crustal thicknesses and 
dynamic compensation mechanisms. 

Implications for Other Regions of Venus 

Other regions of high topography are also candidates for 
gravitational spreading. Even areas that have relatively low 
relief are predicted to undergo gravitational relaxation on 
geologic time scales. Magellan images show evidence of 
extension, which is of modest magnitude and is the latest 
stage of deformation, in other highland and tessera terrains. 

At Alpha Regio, small, widespread graben cut older tectonic 
features, indicating that distributed extension is the youngest 
tectonic event [ deCharon et al., 1991]. Many of the graben 
are parallel to the edges of the elevated terrain of Alpha 
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Regio, suggesting that extension may be a result of 
gravitational spreading [Bindschadler et al. , this issue]. 
Numerous small troughs, interpreted as graben, postdate 
apparently compressional features in Ovda Regio [Saunders 
tl a/., 1991] and could also result from gravitational 
spreading of this equatorial highland. Gravitational 
spreading is also thought to be an important process in the 
evolution of coronae on Venus [Squyres et a/., this issue]. 

Comparison With Other Studies 

In general, the time scales of viscous relaxation predicted 
by these models compare well with those from other studies 
of Venus. For example, in this study, we find that a 14cm- 
high plateau with a surface slope of 3°, a crustal thickness of 
20 km, a thermal gradient of 15 K/km, a basal temperature of 
1040 K, and a model width (plateau, scarp, and plains) of 150 
km relaxes by 25% in approximately 0.1 m.y. In a study of 
viscous relaxation of impact craters on Venus, Grimm and 
Solomon [1988] found that the relief for a 140-km -diameter 
crater would be nearly completely relaxed within 50 m.y., 
given a 15-km-thick diabase mist over an olivine-rich mantle 
and a thermal gradient of 15 K/km. Our results are also in 
general agreement with one-dimensional models of 
gravitational spreading on Venus, which show that viscous 
deformation can occur at geologically significant strain rates 
[ Weertman , 1979; Smrekar and Phillips , 1988]. One would 
expect viscous relaxation on Earth to proceed more slowly 
than on Venus for a given rheology because of the difference 
in surface temperature. However, in models of the relaxation 
of isostatically compensated topography of a given 
wavelength. Bird [1991] found that relief with an amplitude 
of 1 km (2 1cm total relief), a 35-km-thick crust with ductile 
strength given by a diabase flow law, a horizontal 
wavelength of 100 km (surface slope -2°), and a basal 
temperature of 980 K relaxed by 50% within 2000 yr. This 
time scale for relaxation is shorter than would be predicted 
by our models for comparable length scales but higher 
surface temperature, probably because Bird [1991] included 
crustal compensation and had a very different model 
formulation. 

Conclusions 

Magellan radar images display abundant evidence for 
extensional faulting in regions of high elevation ami steep 
topographic slope within Ishtar Terra. We interpret 
extensional faults oriented generally perpendicular to the 
downslope direction to be a result of gravitational spreading. 
Extensional faults that are parallel to the apparent direction 
of shortening in mountain belts are interpreted as evidence 
for gravitational spreading within an actively converging 
orogen. Regions of likely gravitational spreading occur in all 
of the mountain belts surrounding Ishtar Terra and on the 
plateau margins of Lakshmi Planum. Horizontal surface 
strains in these areas are estimated to be 1-20%. 

Volcanism is associated with many of these areas of 
extensional faulting, particularly along the margins of 
Lakshmi Planum and in Danu Montes. Extensional faulting 
cuts older folds and thrust faults in Danu, Akna, Freyja and 
southern Maxwell montes. In Akna, Danu, and Freyja 
montes, some of the extensional faults have been deformed 
by later shortening and possibly horizontal shearing, 
suggesting that extension occurred before the mountain- 
forming process was complete. 


Finite element modeling of the gravitational spreading of 
high topography by flow in the lower crust shows that 
failure and relaxation of relief occur on geologically rapid 
time scales for likely estimates of crustal thickness and 
thermal gradient. The pattern of predicted failure is 
consistent with the observations in Ishtar Terra, in that the 
cumulative horizontal surface strains are low. Extensional 
faulting is predicted to occur first on the plateau, followed by 
extension on the slope, and sometimes by thrusting in the 
plains. Initial failure occurs typically within 0.1 m.y., for a 
thermal gradient of at least 15 K/km. However, an 
accumulated horizontal strain of 1%, which is a more 
plausible threshold for observability in Magellan radar 
images, does not occur typically until 0.01-100 m.y. A 
thermal gradient of 5 or 25 K/km shifts the time scales 
somewhat. 

Viscous relaxation of relief is also predicted to occur 
rapidly. The rates of relaxation depend strongly on the 
temperature (and thus the effective viscosity) at the base of 
the crust, as relaxation occurs primarily through flow in the 
lower crust. For a range of plateau elevations and scarp 
slopes, relaxation of the topographic relief by 25% occurs 
within 10 5 , 10* or 10 10 yr for crustal thicknesses of 30, 20, or 
10 km, respectively, and a thermal gradient of 15 K/km. A 
thermal gradient of 25 K/km gives relaxation times that are 
shorter by up to a factor of 10 3 ; a thermal gradient of 5 K/km 
gives relocation times that are larger by a factor of 10-10 4 . 

These results indicate that the topographic relief in Ishtar 
Terra either is supported by crust considerably stronger than 
indicated by standard flow laws for terrestrial crustal 
materials, or has been built and maintained by tectonic forces 
until times less than millions to tens of millions of years. 
Although an average crater retention age for the high 
northern latitudes on Venus of 50 Ma to 1 Ga [Ivanov et a/., 
1986; Schaber et al. , 1987; Basilevsky et al. t 1987] and an 
average crater retention age of 500 Ma for the entire planet 
[Phillips et al ., 1991a, this issue] tend to favor the 
interpretation that the area is comparatively strong and 
inactive, this age may not be applicable to such small areas as 
mountain belts and plateau margins. Evidence for 
widespread volcanism, which appears to be concurrent with 
extension in many locations, particularly in Lakshmi Planum, 
may be inconsistent with an anomalously cold crust The 
large apparent depth of compensation at Ishtar Terra 
indicates that the long-wavelength topography is probably 
dynamically supported [Grimm and Phillips , 1991] and 
suggests a ready mechanism for recent crustal deformation. 
In the mountain belts, even partial isostatic compensation by 
a crustal root would result in crustal thicknesses greater than 
the largest value used in this study (30 km) and would yield 
an even more rapid relaxation of relief than obtained here. 
We therefore favor the hypothesis that regional deformation 
in the mountain belts and steep plateau margins of Ishtar 
Terra is geologically young (less than 10 m.y.) and may be 
ongoing. 
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^,?a^™H?£ ICAL AND AGE CONSTRAINTS FROM MODELS OF 
GRAVITATIONAL SPREADING IN ISHTAR TERRA, VENUS; Suzanne E. Smrekar 
and Sean C^Solomon, Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts 
Institute of Technology, Cambridge, MA 02144. 


Introducaon. Gravitational spreading is expected to lead to rapid relaxation of high relief due to 
the high surface temperature and associated weak crust on Venus [1,2]. In this study, we use new 
Magellan radar and altimetry data to determine the extent of gravitational relaxation in Ishtar Terra, 
which contains the highest relief on Venus and areas of extremely high topographic slope. Within 
Ishtar Terra, the only mountain belts found on Venus, Akna, Danu, Freyja, and Maxwell, nearly 
encircle the smooth, high (3-4 km) plateau of Lakshmi Planum. Finite element models of this 
process give expected timescales of relaxation and failure at the surface. Using modeling results, 
we attempt to constrain the strength of the crust and timescales of deformation in Ishtar Terra. 
Below we discuss observational evidence for gravitational spreading in Ishtar Terra, results from 
timing 1 of def 1601 “dependent age constraints, and the implications for the rheology and 

Observations . Magellan data have revealed abundant evidence for localized extension 
throughout Ishtar Terra [3,4]. Many of the observed extensional features are oriented 
perpendicular to topographic slope, implying that they formed as a result of gravitational 
spreading. In some of the mountaineous areas, extensional faults occur parallel to the apparent 
direction of shortening. By analogy with active gravitational spreading on Earth [5], we interpret 
this extension to indicate gravitational spreading during convergence. Areas of possible 
gravitational spreading are observed in each of the four mountain belts and along both the southern 
“***“ °f Lakshmi Planum. Horizontal strain in these regions is estimated to be 
1-20%; this value may be an underestimate if blocks are rotated along the normal faults. The 
observational evidence for gravitational spreading in Ishtar Terra is discussed more fully in [4]. 

Models. We use the finite element algorithym TECTON [6] to model the evolution of 
pavitadonal spreading in a vertical section of the crust near die margin of a plateau or the edge of a 
broad mountain belt. Each model includes an elevated plateau 50 km in width, a sloped margin 
aund a lowlands region 80 km in width. The code employs a depth-dependent, viscoelastic 
rheology with non-linear stress dependence and exponential temperature dependence. We adopt a 
Young s modulus of 6x10^ Pa and a Poisson's ratio of 0.25. As a lower bound on the viscous 
strength, a diabase flow law [7]. Although a diabase compostion is most similar to the 
composition determined at Soviet lander sites [8], the flow law determined for diabase may 
underestimate the strength due to partial melting during the experiment [7,9]. For this reason, we 
also employ a websterite flow law [10] as an upper bound on die crustal stmgth. The surface 
temperature is 740 K and increases linearly with depth. Each row of elements in the grid has the 
same viscosity, which is equivalent to assuming that the temperature is constant along both the 
bottom and top of the grid. This results in a somewhat higher thermal gradient in the plains than in 
the plateau. The boundary conditions are zero vertical and horizontal velocity on the bottom of the 
grid, a free surface top boundary, and zero horizontal velocity and free vertical slip at the sides. 

The bottom boundary condition approximates an upper mantle layer that is much stonger than the 
lower crast Details of the models are given in [4].' Brittle failure is evaluated using a Mohr- 
Coulomb criterion. The timing of predicted brittle failure and relaxation of relief are found for 
ranges of plateau height (1-5 km), plateau margin slope (1-30°), crustal thickness (10-30 km), and 
thermal gradient (5-25 K/km). Slopes of 1-30° and plateau or mountain belt heights of 1-6 km are 
observed in Ishtar Terra. Crustal thickness on Venus is predicted to be 10-30 km [1 1 121- the 
average thermal gradient is expected to be 10-25 K/km [13,14]. ’ 

The rate of deformation is largely controlled by the effective viscosity, and thus the temperature 
and stress, at the base of the crust. The greatest deformation occurs at the base of the crust, where 
tiie viscosity is lowest. Temperature at the base of the crust increases with crustal thickness (for a 
fixed thermal gradient), and deviatoric stress near the base of the crust increases with both crustal 
thickness and plateau elevation. The non-linear change in predicted failure and relaxation tim* is 
due both to the exponential dependence on temperature and to the power law dependence on 
deviatoric stress of the effective viscosity. 
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Failure is predicted to occur when the stresses exceed the Mohr-Coulomb criterion. Initial 
failure is located on the plateau. In most cases, failure on the slope occurs shortly after failure on 
the plateau; thrust faulting in the plains occurs typically only for a crustal thickness of 30 km. 
Figure 1 shows the time for the first failure on the plateau as a function of plateau height and 
crustal thickness, for a margin slope of 3° and a thermal gradients of 15 K/km. Models with both 
websterite and diabase flow laws are plotted. No times are shown for models that do not fail 
within 1 Ga. For cases with a diabase rheology, a large plateau height, and a large crustal 
thickness fail almost immediately. A plateau height of 1 km and a crustal thicknesses of 10 km can 
take up to tens of millions of years to fail. With a websterite flow law, no failure occurs unless the 
crust is a least 20 km thick and the plateau elevation is 1-3 km. Changing die thermal gradient by 
10 K/km has a large effect on failure times, with a higher thermal gradient decreasing die failure 
time by a factor of 100-10,000, and a lower thermal gradient increasing times by a similar factor. 
When failure is first predicted, the horizontal surface strain is approximately 0.05%. This strain is 
unlikely to be recognized as normal faulting in Magellan radar images; a strain of 1%, which 
accumulates once significant relxation of the topograpy begins, is probably a more reasonable 
value to compare with observations. As faulting is not explicitly modeled, there is some 
uncertainty in the interpretation of strain after the initial failure. The total strain may be 
underestimated. 

As elevations in Ishtar Terra are among the highest on Venus, we assume that relief has not 
relaxed by more than 25% of its original value. The time at which the average elevation in die 
plateau has relaxed by at least 25% is shown as a function of plateau height and crustal thickness, 
for thermal gradients of 15 and 25 K/km, in Figure 1. For a given thermal gradient and crustal 
thickness, predicted relaxation times for plateau heights of 1, 3, and 5 km vary by a factor of 10 at 
most. A 10 km change in crustal thickness leads to a factor of 10-100 change in relaxation time. 
Larger differences in relaxation time are predicted for models with the same plateau height and 
crustal thickness but with thermal gradients that differ by 10 K/km, as is consistent with the larger 
changes in basal temperature. For a thermal gradient of 25 K/km and a plateau height range of 1-5 
km, all the models are predicted to relax by at least 25% within 100 My for a 10 km thick crust, 

0.1 My for a 20 km thick crust, and 1000 y for a 30 km thick crust A thermal gradient of 15 
K/km gives relaxation times of up to 1 Gy, 100 My, and 0.1 My for crustal thicknesses of 10, 20, 
and 30 km, respectively. 

Discussion. The range of parmeters can be narrowed somewhat by considering the geology and 
morphology of Ishtar Terra. In areas of possible gravitational spreading, the local topographic 
slope ranges from 2-15°. For the puposes of discussion, we will assume a slope of at least 3°. 
Typical local elevation changes are 2-3 km, although elevation changes are as large as 5-6 km in 
Maxwell Montes. Thus plateau elevations as low as 1 km can be eliminated, especially given that 
the topography may have relaxed somewhat from its original value. The very steep topographic 
slopes in many parts of Ishtar Terra suggest that the topography is likely to be at least partially 
isostatically compensated. As mountains in Ishtar Terra rise 4-10 km above the surrounding 
lowlands, the crustal thickness is probably greater than 10 km. Further, the abundance of 
volcanism in Ishtar Terra, as well as the overlap of volcanism and extension in many areas [4], 
implies that the thermal gradient is unlikely to anomalously low. Values of 15 K/km or more are 
probably applicable to Ishtar Terra. Using these limits on the parameter range appropriate for 
Ishtar Terra and assuming a diabase flow law constrains relaxation of the topography by 25% or 
more to occur within 10 My. This result implies that either the crust is anomalous in some way, 
thinner, colder, or stronger than expected, or that mountain building processes have maintained the 
relief in Ishtar Terra until geologically recent times. Next we examine independent information on 
the age of tectonism in Ishtar Terra to attempt to descriminate between these two hypotheses. 

Several types of evidence provide conflicting information on the age of Ishtar Terra. The only 
available measure of absolute age, the average crater retension age for Venus of 500 My, implies 
that much of the surface is quite old [15]. Local geological information must be used to determine 
whether a region is likely to be older or younger than the average age [16]. In Ishtar Terra, die 
abundance of volcansim, especially in Lakshmi Planum, suggests that the surface may be less than 
the average age [17]. In many cases, both the shortening in the mountain belts and areas of 
extension crosscut local volcanism in Lakshmi Planum, implying a younger age for the tectonism 
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[3,4,16]. At least one impact crater in Ishtar Terra appears to be modified by eytenaon due © 
gravitational relaxation [4]. However, the 150 km-diameter crater, Cleopatra, 

Maxwell Montes and appears undefonned. In addition to surface processes, suggest 

that dynamic processed the mantle are actively supporting the long wavelength Wgaptav 
Thescdynanic processes may produce stresses capable of creating the observed short wavelength, 

“teSSKwS Of 4. geologic indicators favor the h^pothras Jb-te 

region is likely to younger than the average crater retension age of Venus, 500 My. On the basis of 
tin? information, aid the assumptions discusses above regarding the appropnatepa^tCT range 
for Ishtar Terra, we intepret the results of our modeling to indicate the tectonic processes have 
maintained the topography until as recently as 10 My, and may still be active. 
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Figure 1 The time at which failure (by normal faulting) first occurs on the plateau is shown as a 
function of plateau elevation for three values of plains crustal thickness. The initial slope of the 
scarp in each model is 3°, and the plains thermal gradient is 15 K/km. Open symbols ^pre^nt 
models with a diabase flow law [X]; filled symbols are for a webstente flow law [X] Times in 
Figures 1 and 2 are determined to the nearest order of magnitude. 

Figure 2. The time at which the height of the plateau has relaxed by 25% is shown as function of 
nlateau height for three values of plains crustal thickness. The lnmal slope of the bounding scarp is 
3°. The open symbols are for a plains thermal gradient of 15 K/km; filled symbols are for a 

thermal gradient of 25 K/km. 
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Radar imaging and altimetiy data from the Magellan mission have revealed a diversity of 
deformational features at a variety of spatial scales on the Venus surface. The plains of Venus 
record a superposition of distinct episodes of deformation, typically manifested in areally 
distributed strain of modest magnitude and coherent over hundreds of kilometers. Ridge belts 
and mountain belts represent successive degrees of lithospheric shortening and crustal 
thickening. As on Earth, mountain belts show widespread evidence for lateral extension and 
collapse both during and following orogeny. Large-scale extensional structures include the 
quasi-circular coronae and broad rises with linear rift zones, both sites of significant volcanic 
flux. Few large-offset strike-slip faults have been observed, but limited local horizontal shear is 
accommodated across many zones of crustal shortening. Several large-scale tectonic features 
have topographic slopes in excess of 20° over a 10-km scale; because of the tendency for such 
slopes to relax by ductile flow in the lower crust, such regions are likely to be tectonically active. 
In general, Venus is tectonically unlike Earth’s oceanic regions, with their large rigid plates and 
deformation confined to narrow plate boundaries, but actively deforming continental regions 
provide a closer analogue. The styles and distribution of tectonic deformation on Venus may be 
consequences of three differences from Earth: the absence of a hydrological cycle and significant 
erosion, strong coupling of mantle convection to the upper mantle portion of the lithosphere, and 
a high surface temperature and thus a significantly shallower onset of ductile behavior in the 


crust. 
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Introduction. While the Pioneer Venus altimeter. Earth-based radar observatories, and the 
Venera 15-16 orbital imaging radars provided views of large-scale tectonic features on Venus at 
ever-increasing resolution [1-3], the radar images from Magellan constitute an improvement in 
resolution of at least an order of magnitude over the best previously available [4], A summary of 
early Magellan observations of tectonic features on Venus has been published [5], but data 
available at that time were restricted to the first month of mapping and represented only about 15% 
of the surface of the planet. Magellan images and altimetry are now available for more than 95% of 
the Venus surface. Thus a more global perspective may be taken on the styles and distribution of 
lithospheric deformation on Venus and their implications for the tectonic history of the planet [6]. 

Generalizations. Tectonic features on Venus are widespread and diverse, and comparatively 
few regions are undeformed. Areas lacking tectonic features at Magellan resolution are restricted to 
relatively young volcanic deposits and the younger impact craters and associated ejecta. Most of 
the surface, during at least some period over approximately the last 500 My [7,8], has experienced 
horizontal strain sufficient to fault or fold near-surface material. 

Tectonic activity on Venus has continued until geologically recent time, and most likely the 
planet is tectonically active at present. Several arguments support this inference. The great relief 
and steep slopes of the mountains and plateau scarps of Ishtar Terra and of the equatorial chasm 
systems are difficult to reconcile with long-term passive support by crustal strength. Because of 
the high surface temperature on Venus, temperatures at which crustal rocks fail by ductile flow 
should be reached at much shallower depths than on Earth [9,10]. Numerical models suggest that 
areas of high relief and steep slope in the Ishtar region should spread under self-gravity by ductile 
flow of the weak lower crust on time scales less than about 10 My [1 1]. Thus the processes that 
build relief and steepen slopes must have been active within the last 10 My. Further, a number of 
features produced by geological processes that have operated more or less steadily during the past 
500 My show evidence of subsequent tectonic activity. About one third of all preserved impact 
craters on Venus have throughgoing faults and fractures, and 1 in 12 are extensively deformed [8]. 
The longest lava channel on the plains of Venus does not progress monotonically to lower 
elevations downstream, indicating that differential vertical motions have occurred since the channel 
was formed [12]. 

Compared with the Earth, horizontal displacements on Venus over the last 500 My have been 
limited. Most of the tectonic features require modest strains and horizontal displacements of no 
more than a few tens to perhaps a few hundreds of kilometers. Plains thousands of kilometers 
across record horizontal strains of order 10~ 2 or less. The great rift systems of Beta and Ada 
Regiones need have extended no more than a few tens of kilometers, on the basis of topographic 
profiles, extended features such as the crater Somerville in Devana Chasma, and analogy with 
continental rifts on Earth [13]. For compressional features, the amount of crustal thickening can be 
estimated from topographic relief and isostatic considerations, but this approach provides only a 
lower bound on horizontal displacements if any crustal material is recycled into the mande at zones 
of underthrusting. For ridge belts 100 km in width and with up to 1 km of relief, horizontal 
displacements of no more than 100 km are required for crustal thicknesses of 10-20 km beneath the 
adjacent plains [14,15]. Mountain belts are exceptional in that greater horizontal displacements are 
required. For a two- to four-fold thickening of the crust beneath the 500-km width of Maxwell 
Montes, the implied minimum horizontal displacement is 1000-2000 km. 

Unlike the Earth, Venus does not show evidence for a global system of nearly rigid plates with 
horizontal dimensions of lOMO 4 km separated by narrow plate boundary zones a few kilometers 
to tens of kilometers across. Predictions prior to Magellan that Aphrodite Terra would show 
features analogous to terrestrial spreading centers and oceanic fracture zones [16] are now seen to 
be incorrect. Evidence for shear is present in the ridge and fracture belts and in the mountain belts, 
but the shearing tends to be broadly distributed and to accompany horizontal stretching or 
shortening. Few clear examples have yet been documented of long, large-offset strike-slip faults 
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J^ E record of impact cratering on 
Venus, recently revealed from high- 
resolubon radar imaging by the Magel- 
lan spacecraft, has been interpreted as 
indicating that the planet underwent 
catastrophic global resurfacing about 500 
million years ago 1 . Although other inter- 
pretabons of the crater characteristics 
have been suggested 2 , the possibility of 
geologically rapid global resurfacing on 
the planet that most resembles our own in 
terms of mass, density and bulk compo- 
sition has generated widespread interest: 
what mechanisms could cause such a 
catastrophe? Several ideas were discus- 
sed at the Fall Meeting of the American 
Geophysical Union last month*. 

The average age of Venus’s surface 
has been estimated from the number of 
larger impact craters (the planet’s thick 
atmosphere has tended to protect it from 
the smaller bodies that would give cra- 
ters smaller than 30 km in diameter) and 
from the cratering rate on Earth and the 
Moon, or from estimates of the number 
of asteroids whose orbits cross that of 
Venus’. The answer, 500 million years 1 - 2 
(a tenth of the age of the Solar System) 
b not in itself remarkable for a planet 
which. Eke the Earth, has a hot, dyna- 
mic ulterior — despite the absence of 
significant erosion. What is remarkable 
is that the craters at all sizes are indis- 
tinguishable from a random population 
(that is, spatial variations in crater 
density are not larger than those due to 
chance) and almost none has been signi- 

firantlir u.. . & 
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volcanic flows beyond the crater rim - 
( *5* even though Magellan images 
show volcanic and tectonic features to be 
nearly ubiquitous on Venus. 

The interpretation of Schaber and 


others' is that most of the surface dates 
from the end of a global resurfacing 
event that ceased about 500 million years 
ago, and that the small fraction of cra- 
tere volcanically embayed or modified by 
deformation indicate that volcanic and 
tectonic activity has since been at much 
lower levels. Phillips and colleagues 2 
have argued that a correlation of loca- 
tions of modified craters with areas of 
low crater density and an inverse cor- 
relation between crater density and 
Magellan radar backscatter (a quantity 
elevated in regions of high topography 
and high roughness, both thought to be 
signatures of comparative geological 
youth) indicate that the Venus surface 
exhibits a spectrum of ages. (The paucity 
of small craters, however, hampers the 
use of crater density in determining the 
relative ages of geological units, a 
method used with the other solid planets 
and satellites.) Phillips et al. have de- 
scribed an alternative model in which 
resurfacing occurs episodically in small 
patches a few hundred kilometres in 
extent, with a characteristic time be- 
tween events of the order of 10 s years 
(Fig- 2). 

Nonetheless, the catastrophic hypoth- 
esis and its possible mechanisms have 
captured the lion’s share of the atten- 
tion. At the meeting, D. L. Turcotte 
(Ojirnen University) proposed that 
global lithospheric overturn operates 
episodicaUy on Venus, and that for the 
past 500 million years the lithosphere has 
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FIG. 1 Magellan radar Image of an impact crater on Venus 
approximately 70 km In diameter. Radar-bright, rough- 
textured ejecta extending up to 60 km from the crater rim 

fhan ^? m ^ rt<able bright flow features extending more 
than 300 tan from the crater walls are thought to date from 
the impact and have not been modified to any significant 
;^l ree ^ later deformation or volcanism outside the crater 
(although the generally radar-dark, smooth floor of the 
crater may contain younger volcanic deposits). 
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Although the mechanism for stabiliza- 
tion of the lithosphere is not an explicit 
component of Turcotte’s scheme he 
argues that a cool and thick lithosphere 

is in better agreement (than one in 

steady-state conductive 
equilibrium) with the large 
values obtained from 
topographic profiles for the 
rigidity of the lithosphere 
at sites of presumed litho- 
spheric bending 4 and with 
the large ratio of 
long-wavelength gravity 
anomaly to long-wave- 
length topographic relief at 
highland areas 5 . 

Parmentier and Hess 6 , 
following Stevenson and 
Bittker 7 , have recently 
suggested that lithospheric 
stabilization may occur on 
Venus because of a de- 
crease in the density of 
lithospheric mantle follow- 
ing melt extraction. But 
they add that the cooling 
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FIG. 2 Two schematic distributions of surface 
ag ® ° n . Venus. With catastrophic 
resurfacing 1 (green histogram), most of the 
surface is about 500 million years old but 
small areas haw younger ages. With episo- 
dic resurfacing 2 (red histogram), there is a 
more Widespread distribution of surface 
ages from 0 to approximately 1.000 million 
years. Both possibilities are broadly 
consistent with the distribution and states 
of preservation of impact craters on Venus 
(Courtesy of R. J. Phillips.) 

of such a layer may subsequently raise its 
density to the point where it becomes 
unstable. Global overturn of this un- 
stable layer, which would sink back into 
the mantle, would be followed by partial 
melting of the upper mantle, global re- 
surfacing, and the gradual development 
of a new buoyant layer and another 
extended period of stability. In simple 
one-dimensional models of this process, 
the time between lithospheric instability 
events is 300-500 million years. 

Others point to time-variable mantle 
convection rather than lithospheric insta- 
bility as a mechanism for global resur- 
facing. An early effort to simulate such 
convection on Venus 8 led to models in 
which the characteristic flow speed and 
mantle heat flux showed large variations 
(by factors of 2-10) at intervals of 100- 

/XX "I 11 ™. y ears - J - Arkani-Hamed 
(McGill University) reported improved 
calculations with better spatial resolution 
which are still time-varying but which 
show significantly weaker fluctuations of 
the order of several tens of a per cent. 

V. Steinbach (Cologne University) drew 
attention to the potential role of upper- 
mantle phase transitions in governing the 
radial character of mantle convection in 
the large terrestrial planets 9 . He and D 
Yuen argue that such phase transitions 
promote the formation of distinct con- 
vening layers in the upper and lower 
mantle, but that as the Rayleigh number 
(governing vigour of flow) in the mantle 
decreases in response to core cooling, 
whole-mantle convection tends to be- 
come favoured over layered convection. 

If Venus cooled more rapidly than the 
arth early in its history, it may now 

Unlon - 
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have a smaller mantle Rayleigh number 
and may have gone through a transition 
to whole-mantle convection, while the 
Earth may still be characterized by 

mantle flow - Such a transition 
would have been accompanied by over- 
ton of the upper mantle, the upward 

Srfodng f hCat ’ and P robab, y b y global 

The lithospheric instability models as 
applied to date for Venus suffer from 
tomg one-dimensional — global syn- 
ouoneity is imposed by assumption, 
tven if such instability mechanisms do 
operate, it is likely that different parts 
of the planet will be at different stages 
in the stabilization and destabilization 
sequence, and whole-planet temporal 
variations will be smoothed out by such 
regional differences. For both litho- 
spheric and mantle-convection mechan- 
isms, attention needs to be given to all 

the terrestrial planets as a group. Both 
the Earth and Mars have gone through 
episodes of great magmatic activity, pre- 
sumably an indication of variable mantle 
convective flux, but neither planet shows 
any evidence of rapid, complete global 
resurfacing during the past 4,000 million 
yeare, and nor do Mercury or the Moon. 

Global resurfacing by large-scale over- 
turn of the lithosphere or upper mantle 
should result in efficient outgassing of 
the upper mantle; consideration of the 
At abundance of Venus’s atmosphere 
suggests that any such widespread out- 
gassing was restricted to times much 
earlier than 500 million years ago 10 
Finally, remember that there are 
alternatives to the catastrophic model 
that seem to be equally consistent with 
the characteristics of impact craters on 
Venus . Further study of the Magellan 
images of impact features and continued 
development of dynamic models of the 
planet s interior should sharpen the com- 
peting hypotheses and may lead to an 
improved understanding of mantle con- 
vection and melting on all of the terrest- 
nal planets, including the Earth. O 

Sean C. Solomon Is In the Department 
of Terrestrial Magnetism. Carnegie In- 
stitution of Washington. Washington DC 
20015. USA. 
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despite evidence from Magelkm images n, 3011 ^ flows external to the crater rim [1 21 

on Venus. One inteipSS features are wid™ ST 

end of a catastrophic global resurfacing event that^eSd m*® SUrface dates fr °m the 

fochon of craters volcanically embayed or modifi<£n^efiZ^ ^ ago ’ *? d Aat * e small 
tectonic activity subsequent to that time that volcanic and 

which resurfacing occurs episodicaTL] S^s a S ^ ****** model, in 

wider spectrum of surface ages akn anrilrc *®tf lew . undretl bolometers m extent and there is a 
crates on Venus (1). A^uX'r of “JS ££ “““‘f 1 wilh ttK characleristics ofS 

?fk e b t en Pmposed, ranging from gSogkX^udd^.n^ f ° f ^^P 1 ^ resurfacing of Venus 
hthosphere [5,6] to strongly destabilization of toe global 

[7-9]. In most of these geophysical models resuSn^l^ gene 1 ra ^^ n m the underlying mande 
yolcamsm. We explore here the SpoffiffSS ?*2? "“P 1 ^ of explicitly by * 

Venus, the primary resurfacing mechanism has been g f eolo S , p alI y recent history of 

We show how such a hypothesis provides at SSn ^ f ° nnadon rather than volcanism. 

observations as do volcanic resurfacing explanation of a wide range of 

tecromc resurfacing hypothesis S? K^^^^“^™P^oas of the 

fnii5 ^ 0b i ervatlons - Any model for resurfacing on VenS^S^^ 1 ^ hlstoi y of the planet, 
following observations: (i) The average crater reinLn™ J?, Id ** consistent with toe 
[1,2,10]. (ii) As noted above toe dSributi^ nf S? ? ona §^? f ? esurface is about5(X)My 
confidence from that of a random population [11 distinguishable at high 

by volcanic flows exterior to toe nm[l 2] Aboit oSh^/*? ° fthe craters ^ embayed 
subsequent to toe impact event [2]; for approxim^lO ? °J Cl ? ers have ^ deformed 
deformation has been extensive n 21 ^ craters the post-impact 

located in areas of low for modi ^d cralere to be 

density and radar backscatter; i.e., smooth dainV^ EL* 5 f Iaverse correlation between crater 
than radar-bright regions of high ODOCTaoh v^n^rfr h ■ l te " de ncy to be more densely cratered 
radar-bright regions on Venus SeTSKtoS com 8 ^ ^ ^ The 
make up large areas of many highland regions^! ll^S P n 5^ed terrain, or tessera, that 

in many lowland plains units [4]. (viii) DeformaHn^nn^r P^^^dy ^ small exposed inliers 
rather than concentrated into naiowSnef^^ SSI v V f t ? ds to ** broad ly distributed 
Platcs[4], (ix) Topography lton^™te^ a ^ oWs y ste ”o ft ec‘onic 
major features have a large gravit^-toroDoeranhvSl^TO^ “i 0 " 8 wavelengths [12]; many 
[13] (x) Evidence for tectonic activity substantial] v mim ^^^PP^cnt depth of compensation 

toe large number of deformed craters includes the <m>»t re ,^ ent , dian 500 M y ago, in addition to 
belte juid plateau scarps of Ishtar Terra and of toe eSTatoriaf/hT^ SteCp sIo Pf s of ^ mountain 
profiles indicating differential vertical movements Sone Sy f?J! J 4, 14 ^ 311(1 elevation 

p e< ? mstaild u PPer mantle of Venus may be stronger toannSw^ivf 18 15 l ^ Evidence that 
Earth and the 450 K greater surface temperature todud^ Sim P le extrapolation from 

craters [2] and large values of elastic lithosphere thickness denvwtf ? ' V U ^ eaxed depths of impact 
topographic profiles across toe margins of coronae ri61 40 ? fl f xural mod els of 

tmosphere is a factor of 4 less than that on Earth as a fraction nf i abun dance of toe Venus 
any widespread outgassing such as might accomoanv law c i!,P fp anet mass t 17 ], suggesting that 
g,6] or upper mantle [8] was restricted to times S lobaI lithosphere 

While episodes of widespread volcanism aUflLSrthL to^S f 5 °° My ag0 [18 ^ & 
have been documented for Mare [19] and Earth r2^n!fninfli! h !? ng ' term average for the planet 
subjected to a global volcanic re/urflctog over Sast 4 Gv ^ pIanets hav ^ ***n 
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of crustal and mantle material are known to vary exponentially with reciprocal temperature so for a 

311(1 t* 1 * * * *® 5 fieId> 111811 rates onflow are expected to be reachedat much 
shallower levels on Venus than on other terrestrial planets. The large values of GTR and apparent 

ilnS ofc °repensation on Venus have been taken as evidence that Venus lacks an asthenosphere or 
pper mantle low- viscosity channel and that mantle convective stresses couple strongly to the 
overlymg lithosphere [21]. In the absence of plate tectonics, these stresses should give rise to 

f^* 155 j r S. adly coherent over large regions. For a sufficiently weak lower 

crust, rates of lower crustal deformation and consequently of surface strain can be high 

* e ® r f f Venus hi 810 /? now preserved, therefore, if the surface temperature was 
comparable to that at present, the higher heat flow associated with early planetary cooling and 
enhanced levels of radiogenic heat production and a mantle convective vigor at least that of the 
o u ed t0 8 e plogically rapid rates of crustal deformation over most, if not all, of 
- SUC 5 T would have characterized by a nearly global extent of complex 
ndged terrain and few impact craters sufficiently undeformed as to be recognizable from surface 
*, some point in the evolution of Venus, however, heat flow will decline to levels 
sufficiently low that the ductile strength of the lower crust will increase rapidly with small 
mcrements of additional cooling. Subsequent to that transition, which might appear to be rapid 
relative to the geological record, rates of deformation will be substantially less, and both volcanic 
deposits and impact craters will persist for long intervals with at most modest deformation of 
iandtorms. The observations enumerated above are consistent with this hypothesis if this 
transition from rapid to modest rates of surface strain accumulation occurred about 500 My ago. 
Implications for Resurfacing History. The tectonic resurfacing hypothesis leads to some 

? at 316 to the resurfacing controversy. If Venus were laterally uniform 

in both crustal thickness and heat flow, then the transition in surface strain rates would occur with 
glob~ synchroneity. That is, there would be a rapid change on a planetary scale from high rates of 
resurfacing to low rates, as is called for by the catastrophic resurfacing model [2], although no true 
catastrophe - and certainly no global outpouring of magma - is involved. While the assumption of 
uniform crustal thickness and heat flow is unreasonable, the unimodal hypsometric distribution for 
Venus suggests that a large fraction of the Venus surface may not depart greatly from this 
assumption; i.e., an apparently "catastrophic" change is not a bad first approximation. Departures 
from a globally uniform change in resurfacing rates are to be expected, however. In particular, 
highland regions, whether they owe their elevations primarily to greater than average crustal 
Sr?S S r^ r t0 f e / ? h j uiced [temperatures at depth, should persist as regions of high strain rate long 
f er ,?f, rat l 0f de £ orm f 1 ? n ^ lowland plains regions has dropped to modest levels. Lowlands 6 
should thus be preferred sites for the preservation of relatively undeformed volcanic deposits and 
impact craters, as is observed [1,3,4]. 

Conclusions. The hypothesis that most resurfacing on Venus has occurred by tectonic rather 
than volcanic processes can account for many of the important characteristics of the planet The 
unusual cratering record on Venus is seen in this light to be a consequence primarily of the 
atmospheric greenhouse and the effect of the high surface temperature on the rheology of the crust 
TtefaxpoAeras leads to the view that the resurfacing history should contain elements of both the 

catastrophic and episodic scenarios for crater removal, with approximately coeval stabilization 

of lithosphere beneath plains regions but more recent tectonic activity concentrated in highlands. 
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The Resurfacing Controversy for Venus: An Overview 

and a Mechanistic Perspective 

Sg3D C. Solomon (DTM, Carnegie Institution of Washington, 5241 

Broad Branch Road, N.W., Washington, DC 20015) 

Two remarkable aspects of the population of impact craters on 
Venus are that craters at all sizes are indistinguishable from a random 
population and that most craters have not been significantly modified 
by tectonic strain or by volcanic flows external to the crater rim, 
despite evidence from Magellan images that volcanic and tectonic 
features are widespread on Venus. One interpretation of these 
observations is that most of the surface dates from the end of a 
catastrophic global resurfacing event that ceased about 500 My ago, 
and that the small fraction of craters volcanically embayed or modified 
by deformation indicate that volcanic and tectonic activity subsequent 
to that time has been at much lower levels. A competing scenario, in 
which resurfacing occurs episodically in patches a few hundred 
kilometers in extent and there is a wider spectrum of surface ages, also 
appears to be consistent with the characteristics of impact craters on 
Venus. While geological and statistical studies of the crater population 
on. Venus offer some promise for distinguishing between these two 
hypotheses, consideration of the possible mechanisms of catastrophic 
and episodic resurfacing provides an independent perspective. 
Potential mechanisms for catastrophic resurfacing of Venus range 
from geologically sudden convective destabilization of the global 
lithosphere to strongly time-dependent heat flux and melt generation in 
the underlying mantle. For most of these mechanisms, resurfacing 
occurs implicitly or explicitly by volcanism. An alternative hypothesis 
is that, at least in the geologically recent history of Venus, the primary 
resurfacing mechanism has been tectonic deformation rather than 
volcanism. Because the rate of surface strain should be controlled by 
the temperature-dependent strength of the lower crust, a geologically 
rapid transition in surface strain rates should be the natural result of 
planetary cooling. This transition would occur at comparable times 
for areas of similar crustal thickness and heat flow (e.g., plains), but 
would be delayed for regions of thicker or hotter crust (highlands). 
The end of the era of high rates of tectonic resurfacing could thus 
appear as a "catastrophe" over the 80% of the planet with an elevation 
within 1 km of the mean, while continued deformation would give rise 
to "episodic" resurfacing to much younger times in the highlands, a 
result consistent with lower crater densities seen in highland regions. 
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The Deformation Belts of Lavinia Planitia 
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Introduction 

Of the many linear tectonic features known to be 
present on Venus prior to the Magellan mission, most 
attention was devoted to the spectacular mountain 
belts found in the Maxwell Montes and elsewhere le.g., 
Campbell et al., 1983; Barsukov et ai, 1986; Bastlevsky, 
19861 It has also been known, however, that there 
is considerable evidence for tectonism on the venusian 
plains. The “ridge belts” described on the basis of 
Arecibo and Venera 15/16 data are particularly good 
examples. In Venera 15/16 images they are most abun- 
dant in Atalanta Planitia, a lowland region in the north- 
ern hemisphere. Like the highland mountain belts, 
ridge belts were inferred to consist of subparallel ridges 
and grooves. However, they appeared to be elevated 
comparatively little above the surrounding plains. The 
ridges and grooves that form them are approximately 
parallel to the margins of the belts, and can be continu- 
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ous along strike for hundreds of kilometers. Spacings of 
individual ridges, as determined from the low-resolution 
Arecibo and Venera 15/16 data, were seen to be typ- 
ically a few tens of kilometers. The belts in Atalanta 
Planitia are themselves crudely parallel, with typical 
belt-to-belt spacings of 300-400 km. In some locations, 
both in Atalanta and elsewhere, individual ridges and 
furrows not organized into belts also were observed. A 
compressions! origin for the ridge belts was the mos 
widely accepted hypothesis prior to Magellan [eg., Bar- 
sukov et ai, 1986; Crump ler et ai, 1986], although 
an extensional origin also was suggested [Sukhanov and 

Pronin, 1989]. , . 

The first deformation belts on the venusian plains 
observed by the Magellan spacecraft have been those 
on Lavinia Planitia, in Venus’s southern hemisphere. 
Solomon et ai [1991] presented a preliminary descrip- 
tion of some of the belts of Lavinia. They noted that 
there were two classes of belts and concluded that both 
had a compressional origin. In this paper we discuss the 
belts of Lavinia in detail, using improved data and ex- 
panded imaging coverage. Because the Magellan images 
have shown the belts to have significantly greater mor- 
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phologic diversity than was apparent in previous data elsewhere presented by Solomon ct ai [1991] are incor- 
sets, we have followed Solomon ct ai (1991] by adopting rect at small spatial scales, and we present improved 
the general term “deformation belts’’ for these features, altimetric data here, 
retaining the term “ridge belt” for one of the morpho- 

logic subtypes. Geologic Units 

In preparing this paper, we have primarily made use 

of three Magellan data products: full-resolution mo- Most of Lavinia Planitia is seen in the Magellan Cl- 
saicked image data records (F-MIDRs), compressed mo- MIDR 45S350, shown in Figure 1. Several geologic 
saicked image data records (Cl-MIDRs), and global to- units, defined on the basis of radar brightness, small- 
pogr&phic data records (GTDRs). F-MIDRs are syn- scale texture, characteristics and abundance of super- 
thetic aperture radar (SAR) image mosaics that retain posed structural features, and apparent relative ages, 
the full spatial resolution of the radar system, sampled can be distinguished in this area. Because the SAR inci- 
at 75 m/pixel. Cl-MIDRs contain the SAR data re- dence angle in this region averages about 25®, the radar 
sampled at 225 m/pixel. Altimetry is taken from GT- return relates primarily to wavelength-scale roughness 
DRs, which contain measurements of local planetary ra- and, less commonly, to meter-scale and larger slopes; 
dius made by the Magellan altimeter, sampled to a uni- the effects of compositional variations should be minor, 
form planetary grid and spatially filtered for smoothing. Consequently, many, and perhaps most, of the geologi- 
The aJong-track spatial resolution of the altimeter varies cal units are distinguished by differences in surface char- 
from about ~8 km near the equator to ~20 km at high acteristics rather than differences in inherent materia] 
latitudes, and the across-track resolution varies simi- properties. For some units, such as lobate flows, radar 
larly from ~I3 to ~30 km [Pcticngill ct a/., 1991]. Due units very likely correspond to those defined using more 
to a software error, the altimetric results for Lavinia and conventional geological criteria. Other units, however, 



Fie. 1. Magellan radar mosaic of Lavinia Planitia, in sinusoidal equal- area projection. The image is centered at 
45® S, 350° E, and the width of the image is 1843 1cm. Letters refer to features discussed in the text. Numbers 
give the locations of other figures. As for all other radar mosaics in this paper, north is at the top and the radar 
illumination direction is from the left. 
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may represent similar materials that have experienced 
different structural evolutions, with consequently differ- 
ent surface textures at centimeter and meter scales. 

We divide the terrains of Lavinia Planitia into four 
broad groups: complex ridged terrains (CRT), textured 
terrains, regional plains, and digitate plains. We have 
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Fig. 2. Complex ridged terrain, exposed through younger plains 
materials. Width of the image is 190 km. See Figure 1 for loce- 


defined one CRT unit, one textured unit, two regional 
plains units, and one digitate plains unit. The follow- 
ing paragraphs describe the units in approximate order 
from oldest to youngest. 

Complex Ridged Terrain 

CRT occurs in Lavinia Planitia as scattered inliers 
surrounded and embayed by younger plains deposits 
(Figure 2). These inliers are characterized by a strong 
linear trend of kilometer-scale or narrower ridges and 
troughs, with a weaker secondary trend of lineaments 
at a high angle to the strong primary trend present lo- 
cally. CRT is much brighter than surrounding plains 
on the SAR images unless the plains are intensely frac- 
tured. 

Textured Terrain 

Many areas in Lavinia are as radar-bright as CRT 
but do not exhibit the same kilometer-scale ridge-and- 
trough pattern. In places, it is possible to resolve a 
very regular fabric of bright lines at a scale of about 
400 m, and it is this fabric that gives this material its 
textured appearance (Figure 3). At the scale of the 
Magellan SAR images, this fabric can be considered 
penetrative. The structural cause of the fabric cannot 
be determined, but a uniformly spaced set of joints is 
one possibility. Textured terrain occurs in two distinct 
physiographic forms: (1) as plains that are significantly 
brighter than adjacent regional plains materials, and (2) 
as long, narrow ridges generally grouped together into 
ridge belts (see below). Crosscutting and embayment 
relationships in some instances indicate that textured 
terrain is younger than CRT and older than regional 
plains. Nevertheless, both textured terrain and regional 



Fig. 3. Textured terrain. Width of the image is 150 km. See Figure 1 for location. 
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plains are of variable brightness, and therefore there are 
places where it can be difficult to define the contact be- 
tween them. The penetrative fabric and, indeed, the 
overall textured aspect are not evident in places where 
textured terrain is intermediate in character between 
well-defined textured terrain and bright mottled plains. 

Regional Plains 

The areally dominant terrain type in Lavinia Planitia 
consists of moderately radar-dark to moderately radar- 
bright plains. Regional plains materials are most likely 
of volcanic origin, but morphologic forms resembling 
flows are rare. 

Mottled plains. These are moderately bright to locally 


bright plains with a “blotchy” or mottled texture at 
kilometer and smaller scales (Figure 4). They are gen- 
erally characterized by abundant superposed wrinkle 
ridges (see below). In places, mottled plains include 
large numbers of small shields or domes from a few to 
about 20 km in diameter. Locally, it is possible to de- 
fine two members of mottled plains that differ in radar 
brightness. At many such locations it is the brighter 
and apparently older of the two members that includes 
the small shields and domes, which are surrounded by 
the darker and apparently younger member (Figure 4). 

Dark plains. These are plains that are darker, have 
more uniform radar badcscatter, and have fewer super- 



Fig. 4. Mottled plains materials. Darker member in the lower right appears to surround small domes and shields 
in the lighter member in the upper left. Also note abundant wrinkle ridges. Width of the image is 155 km. See 
Figure 1 for location. 
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posed wrinkle ridges than mottled plains. Dark plains 
commonly occur as small patches adjacent to or within 
fracture belts (see below), where they are superposed 
on the fractures. In other places, however, dark plains 
appear to grade into the dark member of the mottled 
plains (Figure 4). 

Digitate plains 

This unit consists of complexes of digitate to locally 
lobate flows, many attaining lengths of hundreds of kilo- 
meters (A in Figure 1). Most of these flows are brighter 
on the SAR image, and hence probably rougher at cen- 
timeter scales, than all other plains units. Feeder chan- 
nels (or collapsed lava tubes) are apparent in several 
places, and many flow complexes can be traced ups- 
lope to sources of apparent volcanic origin. Superposed 
wrinkle ridges are sparse, and flow patterns clearly in- 
dicate that adjacent deformation belts already existed 
as topographic ridges when these flows were emplaced. 

Tectonic Features 

Wrinkle Ridges 

Most of the plains areas of Lavinia contain long, nar- 
row, sinuous linear features that appear brighter than 
the background plains (Figure 4). These features are 
generally less than 1 km wide and a few tens of kilome- 
ters long, although some reach widths of several kilome- 
ters and lengths in excess of 100 km. Typical spacings 
range from several to about 20 km and locally up to 


50 km. Where the topography of these features can be 
inferred from brightness variations in the SAR images, 
they are clearly seen to have positive relief, although 
most appear as simple sinuous bright lineaments. We 
interpret these features as compressions! ridges, anal- 
ogous to the wrinkle ridges commonly seen on the lu- 
nar maria and Martian ridged plains [e.g., Plescia and 
Golombek , 1986; Watters, 1988]. 

The wrinkle ridges of Lavinia Planitia show strong 
preferred orientations within subareas of the region 
(Figure 5a). A very strong NNE-SSW trend is apparent 
in the western portion of the map, while the trend is N- 
S in the central portion and NNW-3SE in the northeast 
portion. These trends are also evident in rose diagrams 
of ridge segment orientations (Figure 5 b). 

Wrinkle ridges appear to have formed throughout 
the evolution of plains, because all plains units have 
them, including digitate plains. However, there is a 
clear decrease in the number of ridges from mottled 
plains to dark plains to digitate plains; that is, there 
is a decrease in ridge concentration with decreasing in- 
ferred age. Furthermore, the trends of the ridges ap- 
pear related to locality rather than stratigraphic age 
because the dominant trend of the sparse ridges on dig- 
itate plains is similar to the dominant trend of the more 
abundant ridges on adjacent regional plains. These re- 
lationships suggest that ridge formation is an ongoing 
process. related to local stress fields that did not change 
significantly during the deposition of the various plains 
materials. Although local stress fields varied geograph- 
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Fig. 5b. Rose diagrams of wrinkle ridge orientations in tlx sub- 
areas of Figure 1. Diagrams are normalized to have equal total 
lengths. 


ically during plains formation, they apparently did not 
vary substantially with time in any one locality. 

Grooves 

Over much of the western part of Lavinia Planitia are 
WNW trending radar-bright lineaments. They are geo- 
metrically distinct from the wrinkle ridges in that they 


are both longer and substantially straighter. In many 
instances they are too narrow to be resolved as anything 
but radar-bright lines, but where they can be resolved 
they are seen to be narrow grooves (Figure 6). Lengths 
from 25 to 75 km are common in some areas, while 
lengths of 150 to 200 km are common in others. A few 
grooves have lengths exceeding 500 km. Spacings are 
also variable. The longest grooves (lengths exceeding 
200 km) show typical spacings of 30 to 100 km. Heav- 
ily deformed regions show typical spacings of 1 to 3 km. 
Groove widths are generally about 300 to 700 m. There 
are no instances of these grooves cutting digitate flows, 
but it is not clear if this is because the flows are younger 
or simply because these grooves were never formed in 
the eastern portions of Lavinia where such flows occur. 
The grooves clearly are younger than mottled and dark 
plains. Abutting relationships with wrinkle ridges sug- 
gest that some ridges are older than adjacent grooves, 
some younger. We interpret these grooves to be narrow 
grabens. 

A map of grooves in Lavinia Planitia shows the par- 
allel nature of the grooves in the western part of the 
region (Figure 7). It is noteworthy that in all areas 
where grooves are present, their orientations are per- 
pendicular to those of the wrinkle ridges. Thus, a sig- 
nificant section of Lavinia Planitia is characterized by 
a distinctive pattern of orthogonal compressional and 
extensions! features on the plains between the defor- 
mation belts. 
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Fig. 7a. Map of grooves (lines on plains) and deformation belts in Figure 1. For fracture belts, fracture trends 
are shown by dashed line patterns; for ridge belts, ridge trends are shown by lines with diamond symbols. Shaded 
regions are complex ridged terrain. 



Fig. 7b. Rose diagrams of groove orientations in six subareas of 
Figure 1. Normalization is as in Figure 5b. 


Riigt Belts 

The most obvious tectonic features of Lavinia Plani- 
tia are the prominent deformation belts that transect 
much of the surface. These belts take two very distinct 
forms, which we call ridge belts and fracture belts. An 
example of a ridge belt is shown in Figure 8. The ridges 
that compose them can vary considerably in width both 
from one to the next in a given belt and along the 
length of a single ridge. The maximum width observed 


is about 10 km, and more typical widths are a few kilo- 
meters. Some ridges appear symmetric in cross section, 
while others appear steeper on one flank than the other. 
A smooth, archlike profile is most common, but some 
ridges show a narrow, rugged secondary ridge superim- 
posed on the crest or on one flank of the arch, as is 
common for lunar wrinkle ridges. The ridges are typi- 
cally rather sinuous and commonly bifurcate and merge 
along strike, producing a complex anastomosing pat- 
tern. Some ridges appear to consist of textured terrain 
material, but in most places the unit involved is not 
apparent. 

. It would be desirable to construct detailed topo- 
graphic profiles across ridge belts, but the individual 
ridges are too narrow to be resolved by the Magellan 
altimeter, and the single SAR view obtained during 
the first mapping cycle does not permit stereo analy- 
sis. There is, however, one approach that can be used 
to construct such profiles. One of the belts (Figure 8) is 
cut by several of the tectonic grooves discussed above. 
While these grooves appear very straight where they 
lie on the adjacent plains, they take on a distinctive 
cuspate curvature as they cross each ridge in the belt. 
It is likely that this curvature is an effect of the radar 
viewing geometry; the low radar incidence angle at this 
latitude (about 23 d ) results in the groove, as it passes 
across the elevated ridge, appearing displaced in the di- 
rection of the spacecraft. 

If we make the assumption that the grooves are ac- 
tually straight in map view as they cross the ridges, 
we can use the apparent local deviation of the groove 


13,586 


Squyres bt al.: Plains Tcctonism on Venus 


the image is 205 km. S*« Figure 1 far lotion. A and B ere endpoinU of the 



trend, plus & knowledge of the radar viewing geome- 
try, to construct a topographic profile. Such a profile 
is shown in Figure 9. Most of the ridges in this profile 
are fairly symmetric, with archlike profiles, widths of 
several kilometers, and heights of 200-300 m. (This ap- 
proach to producing a topographic profiles is, of course, 
fundamentally limited by the correctness of the assump- 
tions on which it is based. The detailed heights and 
shapes of the ridges in the profile may be fairly accu- 
rate, since grooves elsewhere tend to be very straight 
on scales comparable to the ridge width. The overall 
height of the belt in the profile, however, should not 
be regarded as accurate, since subtle curvature of the 
grooves over length scales of many tens of kilometers is 
common.) 

The individual ridges within ridge belts differ little 
in appearance from the largest wrinkle ridges, except 
that they tend to be less sinuous. Because both types 
of ridges are inferred to result from compressional de- 
formation, this similarity is not too surprising. Where 
ridge belts deforming textured terrain are embayed by 
plains materials with wrinkle ridges, large, young wrin- 
kle ridges can be distinguished from older ridges of the 
belt. The possibility must be entertained, however, that 
ridge belts can evolve from wrinkle ridges and thus be 
younger than the adjacent regional plains at some locar 
tions. 


Planitia is the fracture belt (Figure 10). In contrast 
to ridge belts, fracture belts are dominated by a com- 
plex pattern of linear to arcuate faults and fractures. 
Many faults appear singly, but others are paired to 
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form grooves with widths ranging from a few kilometers 
down to the resolution limit of the images. The faults 
commonly display complex anastomosing and crossing 
patterns that indicate repeated or progressive deforma- 
tion. Where fracture belts change trend, or where belts 
bifurcate or merge, the pattern of faults is especially 
intricate. 

In general, fracture belts deform mottled plains ma- 
terials but not dark plains materials. The obvious di- 
version of flows and of channels in digitate flow fields 
around some fracture belts in eastern Lavinia (B in Fig- 
ure 1) and the truncation of the fractures by flows in- 
dicate quite clearly that these fracture belts are older 
than the digitate flow fields. 

A noteworthy aspect of a few of the fracture belts is 
that in some parts of them two distinct scales of defor- 
mation are observed. Faults are spaced very closely 
(typically a few hundred meters) almost everywhere 
they are present in fracture belts. However, in a few 
locales these closely spaced faults are concentrated in 
bands of intense deformation that are separated by 
nearly undeformed materials. The bands typically are 
spaced 20-30 km apart. It was reported by Solomon et 
al. [1991] that these bands of intense deformation co- 


incided with local highs in the altimetry, and this was 
interpreted to mean that the faulting might be a conse- 
quence of brittle deformation across the crests of anti- 
clines. However, as mentioned above, this inference was 
based on incorrectly processed altimetric data. In Fig- 
ure 11 we present the same altimetric scan across one 
of the fracture belts, Hippolyta Linea, that was shown 
by Solomon et al [1991], here using the reprocessed 
altimetry. It now appears that the intensely fractured 
bands within the belts are not higher than the adjacent 
materials; if anything, they may be lower. Definitive 
treatment of this problem will have to await acquisition 
and processing of stereographic data, but at present the 
idea that the bands coincide with anticlinal crests ap- 
pears uncertain. An alternative possibility is that the 
bands are regular concentrations of extensional defor- 
mation caused by necking instabilities. 

Much of the faulting that is so prominent in the frac- 
ture belts is similar in appearance, though much greater 
in density, to the faulting on the plains that forms the 
grooves there. In fact, there are a number of instances 
(e.g., Figures 1 and 10 and Figure 17 below) where fault- 
ing in fracture belts continues onto the plains, at much 
lower density, to form widely spaced grooves. As is the 



Fig. 10. A fracture belt. Width of the image is 260 km. See Figure 1 for location. The altimeter ground track far 
orbit 415 is shown. Arrows correspond in position to arrows in Figure 11. 
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Fig. 11. Altimetry d&t& from orbit 415. Arrows correspond m 
position to arrows in Figure 10. 


case for the ridge belts and the wrinkle ridges, then, 
there appears to be a close relationship between fault- 
ing in fracture belts and the formation of grooves. 


Altimetry 

Magellan altimetry data have also provided the first 
detailed look at the topography of Lavinia Planitia. 
Figure 12 presents a topographic contour map for the 
Lavinia region. Overall, Lavinia is a broad, low plain, 
lying at a typical planetary radius of about 6050.5 km. 
The regional topography rises gradually by more than 
1 km to the northwest, northeast, and southeast. The 
general rise in elevation to the east continues well be- 
yond the edge of our map and culminates in a line of 
several large coronae south of Alpha Regio [Stofan el aL, 
this issue]. A prominent characteristic of the altimetry 
of Lavinia is that both types of deformation belts are 
elevated above the surrounding plains. 

The fracture belts appear in general to reach higher 
elevations than the ridge belts. A number of the frac- 
ture belts reach elevations that are as much as 1 km 
above the level of the adjacent plains, and elevations 
500 m above the plains within fracture belts are the 
norm. Particularly high elevations commonly are found 
in complex and intensely deformed regions present at 
bends and junctions in fracture belts. In contrast, max- 
imum relief of 200-300 m is more typical for the ridge 
belts, and some of the ridge belts are not evident at 
all in the altimetric data. Some care must be exer- 
cised in interpretation of the altimetry, however, as de- 
termination of surface elevation for a given altimeter 
pulse can depend sensitively on both the large-scale and 
small-scale surface roughness within the pulse footprint. 
Roughness on a large scale can cause multiple echoes 
that severely complicate the range-finding algorithm. 
Roughness on a small scale can result in a surface scat- 
tering law that deviates from the one assumed in fitting 
the shape of the leading edge of an echo, again affecting 
the range-finding process. In any case where two ter- 
rains of significantly different roughness are observed, 
systematic variations in elevation between them must 
be treated with caution. 



Lavinia: Contour Interval « 100m 
Latitude: 44.ZS - 50.0S 
Longitude: 347.8E - 356.9E (at north edge) 

Fig. 12. Topographic contour map for the area covered by Figure 
1, In sinusoidal equal-area projection. . Elevations are referenced 
to a planetary radius of 6050.0 km. 


Geologic Relationships 

Evidence for Horizontal Shear in Deformation Belts 
Solomon et al. [1991] presented two pieces of evi- 
dence for modest amounts of horizontal shear parallel to 
the axes of deformation belts. These are shown in Fig- 
ures 8 and 10. In Figure 8, an old set of plains grooves 
exhibits an S-shaped bend consistent with distributed 
left-lateral shear as it crosses the ridge belt. In Fig- 
ure 10, the fracture belt exhibits a number of instances 
where paired faults form rhombohedral downdrops, in 
a manner similar to what is seen in terrestrial settings 
where extension is coupled with shear. In this instance, 
the inferred sense of horizontal shear across the belt is 
right-lateral. , 

There are several other instances in Lavinia Plani- 
tia where patterns developed within deformation belts 
can be interpreted as very large-scale analogs of famil- 
iar kinematic indicators seen in terrestrial shear zones. 
One of these (Figure 13) shows a series of S-shaped mi- 
nor folds along the lengths of individual ridges within 
a ridge belt. S folds of this geometry are indicative 
of left-lateral (sinistral) movement and in this instance 
would imply sinistral shearing parallel to the length of 
the ridge belt. The belt showing the S folds termi- 
nates to the SSW against a contact with mottled plains, 
and there is no indication of sinistral slip within these 
plains. Consequently, the inferred belt-parallel sinistral 
slip must be older than the mottled plains. It either rep- 
resents a late stage in the development of the ridge belt 
or implies that the ridge belt is transpressional (com- 
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pression plus sinistra] shear) rather than purely com- 
pressional in origin. 

A second example of a possible kinematic indica- 
tor (Figure 14) involves the distortion of well-developed 
penetrative fabric in textured terrain as the fabric im- 
pinges on a throughgoing, strongly lineated boundary of 
* ridge belt. This boundary abrutly truncates all struc- 
tures within the belt and is interpreted to be a shear 
rone. If this inference is correct, the distortion of the 
fabric at this location implies a right-lateral (dextral) 
sense of shear. 

Age Relationships Among Belts 

The deformation belts in Lavinia did not all form 
contemporaneously. Instead, they show some signifi- 
cant variability in age, both with respect to the plains 
materials surrounding them and with respect to one an- 
other. As noted above, a number of belts, particularly 
in western Lavinia, appear to deform the same mate- 
rial that constitutes the adjacent plains, with little or 
no subsequent volcanism. In the ridge belts, the tran- 
sition from dark or mottled plains adjacent to the belt 
to the textured terrain of the belt is not necessarily a 
stratigraphic one in all cases; we cannot rule out the 
possibility that in some instances the texturing simply 
resulted from the folding of the plains material. How- 
ever, there are other very clear instances where belts are 
embayed and partially buried by lavas that postdate the 
belt deformation. The N-S fracture belt at C in Figure 
1 is one of these, and it may have once been contiguous 
with another aligned and partially buried N-S belt at 

<y. 

An example of a clear crosscutting relationship in- 
volving two belts is shown at D in Figure 1. Here a N-S 
belt has been largely buried by lavas, and these lavas 
have in turn been deformed to form an E-W belt that 
crosscuts the N-S one. This case in particular demon- 
strates that belt formation in Lavinia took place over an 
extended period of time, during which which belt for- 
mation episodes were interleaved with plains volcanic 
events. 

Regional Tectonic Trends 

As noted above, there are two classes of deformation 


Fig. 14. Penetrative fabric in textured terrain that ia distorted 
with a geometry consistent with right- liters] shear across a ridge 
belt. Width of the image is 125 km. See Figure 1 for location. 


belts in the Lavinia region, and each appears related 
to tectonic features on the intervening plains. Specifi- 
city. ridges of the ridge belts are morphologically sim- 
ilar to and in some cases contiguous with plains wrin- 
kle ridges, and faults of the fracture belts are similar 
to and in some cases contiguous with plains grooves. 
In some parts of Lavinia, however, these relationships 
may be more fundamental, as suggested by Figures 5 
and 7. The morphologic form taken by each deforma- 
tion belt, ridge belt or fracture belt, appears related 
in many instances to the orientation of the belt with 
respect to the orthogonal tectonic features on the ad- 
jacent plains. Belts whose orientations are closest to 
those of the plains grooves in their vicinity tend to be 
fracture belts, while those whose orientations are closest 
to those of local wrinkle ridges tend to be ridge belts. 
Even more notably, there are instances, as in Figure 15, 
where a single belt changes from one trend to the other, 
and in so doing also changes tectonic style from one belt 
type to the other. In addition, there are some belts (e.g., 
Figure 16) where both ridges and fractures are common 
within the belt and where the ridges lie parallel to local 
plains wrinkle ridges and the fractures lie orthogonal to 
them and parallel to local plains grooves. 

There are exceptions to this relationship. The most 
obvious one is the N-S fracture belt at C in Figure 1. 
In this same belt, some ridges also are present and lie 
parallel to fractures, rather than perpendicular to them. 
Other belts departing from the general relationship oc- 
cur farther to the east. Despite these exceptions, the 
relationship described above appears sufficiently preva- 
lent, particularly in western Lavinia, that it must be 
accounted for in some fashion. The belts in central and 
eastern Lavinia that deviate from the general pattern 
are in mo6t cases very clearly buried and embayed by 
younger volcanic deposits and perhaps are formed in a 
stress field distinct from the one responsible for most of 
the wrinkle ridges and grooves. 
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Fig. 15. A deformation belt that changes from a ridge belt to a fracture belt along iU length. Width of the image 
is 290 km. See Figure 1 for location. 


Discussion 

Geologic Evolution of Lavinia Planitia 

The description above of geologic units and tectonic 
features in Lavinia Planitia leads us to some preliminary 
interpretations: 

Regional tectonic pattern on the plains . The orthog- 
onal pattern of compressions] wrinkle ridges and exten- 
sions! grooves that transects the plains units over much 
of Lavinia Planitia is remarkable in its regularity and in 
both its spatial and temporal extent. These character- 
istics indicate that at least the near-surface materials 
of Lavinia have been subjected to a stress field that has 
been generally uniform, both spatially and temporally, 
for much of the region’s history. This conclusion is not 
unique to Lavinia, since a variety of regular tectonic 
patterns are observed on the plains elsewhere on the 
planet [Solomon et ai t 1991, this issue]. 

It is tempting to think of the wrinkle ridges and 
grooves on the plains as the simple manifestations of 
thrust faulting and normal faulting, respectively, and 
this interpretation may well be correct. If so, it has 
some interesting implications. A common approach to 
predicting faulting style from stress field descriptions 
is that of Anderson [1951]. With the Anderson fault- 
ing criteria, however, no single set of values of the three 
principal stresses 02 , and <r 3 simultaneously predicts 
orthogonal grabens and thrust faults. Applying An- 
derson criteria rigorously, one can obtain either thrust 


faults or grabens from a given stress field orientation, 
but not both. 

Orthogonal thrust faults and grabens may both be 
formed, according to Anderson criteria, if the two stress 
components in the plane of the surface maintain their 
relative values but vary appropriately with respect to 
the vertical stress. One way this could take place would 
be for the in-plane stresses to vary relative to the ver- 
tical stress over time. Another would be if fault initia- 
tion for the two classes of features took place at differ- 
ent depths [Golomlek t 1985], near the surface for thrust 
faults and at depth for grabens, so that in the former 
case the vertical stress is <r 3 and in the latter case it is 
<r\. (The depth is limited, however, since many grabens 
are very narrow.) 

In either of these scenarios there is a stress state, 
intermediate in either time or depth, for which Ander- 
son criteria would predict strike-slip faulting. Clear ev- 
idence for strike-slip faulting is rare on Venus (and in- 
deed on most of the terrestrial planets save Earth), and 
is lacking entirely on the plains of Lavinia. One expla- 
nation offered for the absence of strike-slip motion on 
other planets is that layered mechanical discontinuities 
concentrate stress, and hence fault initiation, at certain 
depths (and hence at certain confining stresses) where 
faults with styles other than strike-slip tend to be ini- 
tiated [ Golombek , 1985]. Another possibility is simply 
physical restriction of lateral movement. A planet has 
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► A dcfom^ion bell that exhibit* both ridges and fractures, with approximately orthogonal orientations. 
Width of the image is 240 km. See Figure 1 for location. 


a free surface and a ductile region at depth, so thrust 
and normal faulting within the strong upper crust can 
always be accommodated by vertical motions, and the 
stress at the base of a fault cutting the entire strong 
upper layer can relax by ductile flow. On a one-plate 
planet or within the interior of a single lithospheric 
plate, however, there may be restrictions on the extent 
of horizontal strain, and stress relaxation at the ends of 
strike-slip faults will not be comparable to that for dip- 
slip faults. It may be that such restrictions inhibit most 
large-scale (i.e., easily observed) strike-slip motions in 
Lavinia. Schultz and Zuber [1991] have recently sug- 
gested that in situations where the lithosphere is ver- 
tically loaded, tensile failure rather than shear failure 
may relieve much of the stress and may suppress the 
zone in which shear failure criteria predict strike-slip 
faulting. 


Deformation belts . In contrast to pre- Magellan ex- 
pectations, we now know that there are two distinct 
classes of deformation belts (ridge belts and fracture 
belts) in Lavinia. A particularly important question is 
whether or the not ridge belts and fracture belts both 
owe their origin to a similar mechanism. Evidence sug- 
gesting that they do includes geometric similarities be- 
tween the two classes of belts (widths, positive topo- 
graphic signatures, patterns traced across the plains) 
and, particularly, the observation that a few belts un- 
dergo a transformation from ridge morphology to frac- 
ture morphology along their length. Similar changes in 
morphology along the length of a belt occur in the de- 
formation belts of Atalanta Planitia as well [Solomon et 
al. t this issue]. It is clear, however, that whatever the 
underlying cause of belt formation, the surface tectonic 
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expressions of ridge and fracture belts are markedly dif- 
ferent. 

Because of the gross similarities of both types of 
belts, and particularly because both are elevated, one 
reasonable working hypothesis is that both result from 
belt-normal crustal shortening and thickening [Solomon 
et «/., 1991], In the case of the ridge belts, this argu- 
ment is straightforward: the ridges are interpreted to be 
folds formed by ridge-normal compression. In the case 
of the fracture belts, however, the argument is really 
based only on the altimetry. All of the belts lie higher 
than the surrounding plains and hence are attributed 
in this interpretation to a thickening, via belt-normal 
shortening, of the buoyant crust. 

A difficulty with the idea that both types of belts 
share a common origin is that the most straightforward 
interpretation of the observed tectonic features in the 
fracture belts runs directly counter to this interprets^ 
tion. While near-surface materials in the ridge belts 
have been deformed to produce folds, the same mate- 
rials in fracture belts are intensely faulted in what ap- 
pears to be a manifestation of near-surface extension. 
One could interpret these belts, then, as being wholly 
extensional features, with the uplift a consequence of 
lithospheric heating or igneous intrusion. Alternatively, 
in the compress ion al interpretation, the observed fault- 
ing could be the consequence of stretching of brittle 
surface materials across the crest of an archlike uplift 
whose underlying origin is compressions!. The strongest 
argument that both belt types owe their origin to the 
same process may be that single belts exist that exhibit 
both morphologies, depending on their local trend. 

While the great extent and uniformity of the orthog- 
onal tectonic pattern of Lavinia indicate that the near- 
surface materials there have been subjected to a stress 
field of significant spatial and temporal regularity, the 
deformation belts of Lavinia, in contrast, are geometri- 
cally very irregular and can change trend substantially 
along their length, bifurcating and merging to form an 
intricate pattern. If there is a single underlying belt 
formation process, then the irregular orientations of 
the belts with respect to the orthogonal pattern of the 
plains suggest that this process is not related in a sim- 
ple manner to the one responsible for the near-surface 
stresses expressed on the plains. The morphologic mani- 
festation of each belt nonetheless commonly depends on 
the orientation of the belt with respect to the tectonic 
pattern of the plains: ridge belts form where the belt 
is perpendicular to the regional compressive stress, and 
fracture belts form where the belt is perpendicular to 
the regional extensional stress. In both cases there ap- 
parently can be some distributed shear across the belt, 
as indeed is likely to be the case when the belt axis is 
not perfectly orthogonal to any principal stress. 

The source of stress responsible for formation of the 
deformation belts of Lavinia is not known, but man- 
tle convection must be considered the prime candidate. 

It has been suggested that the Venusian mantle may 
lack an Earth-like asthenosphe re [e.g., Phillips , 1990]. If 
this is the case, convective stresses can couple strongly 
to the lithosphere and can lead to deformation of the 
lithosphere at wavelengths comparable to the spacing 
of the Lavinia deformation belts [ Zubtr , 1987, 1990]. 


The surface manifestation of convective stresses, how- 
ever, will also depend on the rheologic structure of the 
crust. A model in common use for crustal and upper 
mantle rheology on Venus involves a crust of basaltic 
composition overlying an olivine-rich mantle. This con- 
figuration can lead to a layered rheologic structure, with 
a strong upper crust, a ductile lower crust, a strong up- 
per mantle, and ductile mantle material below [ Zubtr , 
1987; Bantrit and Golombek , 1988]. A weak lower crust 
may be relatively ineffective in transmitting convective 
stresses to the upper crust. 

We cannot presently answer the question of whether 
the mantle motions and resultant stresses that led to 
formation of the two classes of belts in the Lavinia re- 
gion were fundamentally the same. In some respects 
the interpretation of the data is more straightforward if 
they were not, since the surface morphologic expressions 
of the two belt types are quite different. However, we 
have already noted the arguments for the view that both 
belt types share a common underlying origin, and these 
arguments appear fairly persuasive. Because the pat- 
tern of the belts is so complex, stresses produced in the 
strong upper mantle by convection must also have been 
complex, both spatially and temporally. The stresses 
that produced surface deformation on the plains, how- 
ever, must have been spatially and temporally quite 
regular. The surface stress field must therefore have 
resulted from the superposition of (1) stresses gener- 
ated by mantle convection and transmitted across the 
ductile lower crust with (2) stresses in the strong upper 
crust of unknown origin that were much more regular 
in space and time. The former stresses are expressed at 
the surface primarily by the elevated belts themselves. 
The latter stresses must have dominated the in-plane 
components of stress in near-surface materials, control- 
ling both the orientation of the plains tectonic features 
and the style of deformation within the belts. A mixed 
stress field of this sort might be enabled by a ductile 
lower crust that incompletely transmits stresses from 
the upper mantle to the upper crust. 

Our observations of the Lavinia region raise a number 
of questions: 

1. What process is responsible for the stress field 
that produced the orthogonal tectonic pattern on the 
plains, both here and elsewhere on Venus? 

2. What process is responsible for the stresses that 
produced the elevated deformation belts, and for the 
concentration of the deformation into discrete belts? 

3. What do the dominant wavelengths of deforms 
tion, both within and among the belts, indicate about 
the mechanical properties of the venusian crust and 
mantle? 

4. If the two classes of belts share a common origin, 
is there a fully self-consistent stress distribution and 
history that can produce the deformation observed? 

5. If the two classes of belts have fundamentally dif- 
ferent origins, why are they so similar in some respects? 
Why are there belts with both characteristics, and espe- 
cially belts that change style along their length? Why 
are the fracture belts different from rift zones elsewhere 
on the planet? 

In the two sections that follow, we make a first at- 
tempt to address questions 2 and 3 above. The others 
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will be left for consideration in future work. 
Dynamical Models for Belt Formation 

We have argued above that the large-scale morphol- 
ogy of deformation belts, and of the ridge belts in par- 
ticular, may be consistent with the hypothesis that 
th«e features are products of lithospheric shortening 
and crustal thickening. This was also the conclusion 
reached by most investigators prior to the Magellan mis- 
sion, and previous workers have postulated this shorten- 
ing to be the result of convective downwelling centered 
beneath lowlands regions [Zuber, 1990; Bindschadlcr ci 
al., 1990; Phillips et al., 1991]. In this section we use 
finite element models to investigate quantitatively the 
relationship between mantle convection and crustal de- 
formation, with particular emphasis on the conditions 
favoring the formation of deformation belts. 

The geometry and characteristic spacing of the de- 
formation belts of Atalanta Planitia suggest that these 
features result from lithospheric instabilities induced 
by horizontal extensions! or compressional stress [Za- 
ber, 1987]. There are two characteristic scales of de- 
formation associated with these belts: one describing 
the spacing of features (e.g., ridges) within the belts, 
and one describing the spacing of the belts themselves. 


Linear stability analysis indicates that two scales of de- 
formation can arise given sufficient compression or ex- 
tension if the lithosphere is Theologically stratified, i.e., 
if a weak lower crustal layer is sandwiched between a rel- 
atively strong upper mantle and a strong upper crust. 
A stress-dependent rheology is required for the instabil- 
ities to grow. Znber [1990] suggests that the congres- 
sional stresses required to induce the instabilities may 
be supplied by mantle convection. This hypothesis re- 
quires that the stresses associated with the motions of 
the mantle couple into the overlying crust. The high 
correlation of long-wavelength topography and gravity 
and the large apparent depths of compensation suggest 
that Venus lacks an upper-mantle low-viscosity zone 
[Kiefer et aL, 1986; Bills et al., 1987; Phillips et al, 
1991; Smrekar and Phillips, 1991], If so, then mantle 
convection should indeed couple strongly to the litho- 
sphere, and patterns of mantle flow should have recog- 
nizable signatures in surface topography and deforma- 
tion [Phillips, 1986, 1990]. 

On the basis of analytical models of convection- 
induced crustal flow [e.g., Bindschadler and Parmen- 
tier, 1990; Schmeling and Marquart, 1990], the crust 
may be regarded as having two responses to convective 
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downwelling. The first is a “flexural” response, dur- 
ing which the surface and crust-mantle boundary 
bo thought of as deforming “in phase.” Following this 
mode of deformation, the crust responds by thi cke ning, 
and eventually the surface and the crust-mantle bound- 
ary become anticorrelated. 

Simons et at [1991] have proposed that deforma- 
tion belts like the ones observed in Lavinia Planitia 
can form during the transition from the initial flexural 
response to the later phase of crustal thickening, and 
we expand on this suggestion here. We illustrate this 
proposal by means of idealized finite-element models of 
the dynamic interaction between the crust and man- 
tle. The crust and mantle in these models are treated 
as viscous fluids, and a modified version of the two- 
dimensional Cartesian finite-element program ConMan 
[King et al., 1990] is employed. A passive marker chain 
tracks the crust-mantle interface and permits variation 
in the crustal buoyancy as well as specification of dif- 
ferent crustal and mantle rheologies. In addition to the 
flow field, the stress field in the lithosphere and the sur- 
face topography are readily calculated. 

We use an irregular finite element mesh that is 28 
elements high and 80 elements wide (Figure 17). Our 
maximum resolution is in the 40-km-thick top layer, 
where each element is 4 km high and 10 km wide. In 
all, the mesh is 800 km in the horizontal dimension and 
400 km in the vertical dimension. We impose free-slip 
boundary conditions on the top and side walls, with 
no flow through these walls. On the bottom boundary, 
flow is constrained to be vertical, with no deviatoric 
normal stress permitted. In effect, this last boundary 
condition gives us a virtual 800 km by 800 km box. 
The surface topography is calculated from the vertical 
stresses on the top wall of the box. Top and bottom 
temperatures are fixed at 500 # C and 1250°C, respec- 
tively. Initially, we impose a linear temperature gradi- 
ent across the thermal boundary layer and set the rest 
of the mantle to be nearly isothermal with a small si- 
nusoidal temperature perturbation imposed throughout 
the box to initiate flow. 

Our exploratory models have a constant viscosity 
mantle overlain by a crust of greater constant viscosity. 
As noted above, the viscosity in the actual Venusian 
crust and mantle is likely to vary with depth; in these 
simple models the higher viscosity crust is intended to 
represent a strong surface layer capable of being de- 
formed by stresses induced by mantle convective mo- 
tions. We adopt a reference mantle viscosity of 10 21 
Pa s, a crustal density of 3000 kg m -3 , and a mantle 
density of 3300 kg m -3 . We use initial crustal thick- 
nesses of 12 and 25 km, and ratios of crustal viscosity 
to mantle viscosity of 25, 50, and 100. In all these mod- 
els, convection produces horizontal compressive stresses 
in the lithosphere that overlies downwelling mantle and 
extensions] stresses in the lithosphere that overlies up- 
welling mantle. As the convective vigor increases, so 
does the magnitude of the compressive stress over the 
downwelling, with stress in the crust reaching values 
in excess of 100 MPa in about 100 m.y. (Figure 18). 
Note that the rate of increase in compressive stress de- 
creases with increasing crustal viscosity. This is because 
the stronger the crust, the more the development of the 
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Fig. 18. Development of horizontal comprestive stress versus time 
in models of convectively induced crustal deformation. Stress is 
measured at the top right-hand (downwelling) side of the model. 
Models shown include crust/mantle viscosity ratios of 25 (solid 
line), 50 (dashed Une), and 100 (dotted line), and a crustal thick- 
ness H of 12 or 25 1cm. Note the similarities among the curves 
and the peak in oompressional stress. 


convective instability in the mantle driving the deforma- 
tion is impeded. We also find that the magnitude of the 
peak compressive stress achieved above the downwelling 
increases with higher viscosities and/or with thinner ini- 
tial crustal layers; the stronger the crustal lid, the larger 
the tractions from mantle convection supported by the 
crust. Since force balance on the crust requires that 
shear traction integrated along the base be balanced 
by normal stresses integrated through its thickness, the 
thinner the crust, the larger the horizontal stresses. 

The topography over the downwelling undergoes a 
complex evolution in these models (Figure 19). Initially, 
the effect of stresses induced by downgoing mantle dom- 
inates, producing topographic lows. After reaching a 
minimum, the topography grows due to the isostatic 
contribution from the thickened crust. We suggest that 
if the deformation belts are the result of compressive 
stress and that mantle convection is the source of this 
stress, then the belts will develop at or before the time of 
maximum horizontal compressive stress. This time cor- 
responds to the transition from decreasing to increasing 
topography. 

Theory derived by Biot [1961] shows that horizontal 
compression of a buoyant viscous plate overlying a vis- 
cous substratum will have a dominant wavelength that 
depends on the difference in layer densities A p, the ratio 
of viscosities the applied stress P t and the layer 

thickness H c . For values of the quantity P/(ApgH c ) 
less than about 7, the amplification of folds becomes 
insignificant. For one of our models described above, 
with fjc/Tm = 100, H c = 12 km, A p s 300 kg/m 3 , and 
P = 100 MPa, the dominant wavelength is about 100 
km. However, these parameters are near the threshold 
conditions for significant fold amplitudes. 

Our initial results show that coupling of convective 
motions to the lithosphere will cause significant crustal 
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Fig. 19. Evolution of surface elevation (top) and the depth to the crust-mantle interface (bottom) with time for 
* model ofcnistal deformation induced by mantle convection. Note the initial decrease in elevation over the tone 
of downwelling (downward arrow), followed by an eventual increase in elevation. Note also the steady increase in 
cnistsJ thickness over the tone of downwelling. The initial crustal thickness, H c , is 12 km, and the ratio of crustal 
viscosity to mantle viscosity is 70. 


deformation. This deformation is complex and highly 
time dependent. Our results are preliminary, and mod- 
eling with more realistic rheologies is in progress. How- 
ever, the simple rheologies assumed here are sufficient 
to generate stress fields that, together with simple the- 
ory [Biot, 1961], predict the predominant wavelength of 
deformation seen in Lavinia Planitia. This theory may 
also explain why the belts are only generated for a fi- 
nite period of time before they become overwhelmed by 
larger-scale deformation. 

Buckling Within Individual Ridge Belts 

We have argued that the ridge belts observed on 
Lavinia Planitia are fold belts, formed by buckling of a 
strong upper crust in response to ridge-normal compres- 
sion. If this interpretation is correct, then we may use 
observations of ridge geometry, along with a simple me- 
chanical model of buckling, to place some quantitative 


constraints on the structure of Venus’s crust at the time 
of ridge formation. We model the ridges here as having 
formed due to a buckling instability in an elastic layer 
with a depth-dependent yield strength subjected to a 
compressive end force. Our genera] approach follows 
that of Banerit and Golombek [1988]. We consider the 
strength of the venusian lithosphere to be determined by 
two deformation mechanisms: brittle failure and ductile 
flow [see Brace and Kohlstedt , 1980]. At a given depth, 
the local pressure and temperature determine the yield 
stress for each mechanism; the lithospheric strength is 
taken to be the weaker of the two. 

In our modeling, the critical differential stress for 
brittle failure is taken to be r & /i<r n [Byerlee, 1978], 
where r and a n are the shear and normal stresses, re- 
spectively, and fi is the coefficient of friction. The criti- 
cal differential stress for ductile flow is derived from the 
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creep equation i = A<r n exp(— Q/PT), where c is the 
strain rate, a is the differential stress, R is the gas con- 
stant, T is the absolute temperature, and Q and n are 
experimentally determined constants. We take values 
for these constants appropriate for dry diabase [Carxs- 
tan, 1982], 

An example of a lithospheric strength profile for hor- 
izontal compression calculated from the above models is 
given by the heavy solid line in Figure 20. The strength 
increases linearly with depth near the surface (in the 
brittle regime) due to the increasing lithostatic stress 
and then decreases exponentially (in the ductile regime) 
due to the thermal gradient. For this figure we use a sur- 
face temperature of 720 K, a thermal gradient (assumed 
constant) of 15 K km -1 , and a strain rate of 10“ 15 s* 1 . 
Now, consider how the strength profile of Figure 20 re- 
sponds to an applied stress of of 70 MPa (the dashed 
vertical line in the figure). Only those depths with a 
lithospheric strength exceeding the applied stress (be- 
tween Z x and Zb) will support the entire stress load 
elastically. Depths above Z t will only support part of 
the load elastically (up to the yield stress) and will un- 
dergo brittle failure. Depths below Z\ will only support 
part of the load and will undergo ductile creep. 

At the crust-mantle boundary, a change in compo- 
sition can create another rheologic transition. Mantle 
material on Venus has been assumed in most models 
to be olivine-rich in composition and hence to possess 
much higher yield stresses than the crust, leading to an- 
other brittle-ductile failure envelope in the upper man- 
tle [e.g., Zuber , 1987; Banerdt and Golombek , 1988]. In 
line with the calculations above, this deeper strong zone 
could be responsible for longer wavelengths of deforma- 
tion on Venus. Here we only consider short-wavelength 
folding and take it to be a function only of the crustal 



Fig. 20. A typical strength envelope for the venusian crust, show- 
ing strength as a function of depth, with an upper regime con- 
trolled by frictional failure on faults and a lower regime controlled 
by ductile flow. Calculations are for a dry diabase rheology, a sur- 
face temperature of 720 K , a thermal gradient of 15 K km*" 1 , and 
a strain rate of 10“ 15 s* 1 . 


failure envelope, neglecting the potential influence that 
the coupling between the crust and mantle discussed in 
the previous section could have on the stress field. 

An elastic plate overlying an inviscid fluid will buckle 
if subjected to a sufficiently large compressive end load. 
The critical force required is given by [Turtotte and 
Schubert, 1982, p. 124] 



where E t h , and v are the Young’s modulus, thickness, 
and Poisson’s ratio of the plate, respectively, p is the 
fluid density, and g is the gravitational acceleration. 
Given this force, the elastic plate will deform with a 
characteristic wavelength 

IwSd" » 

This simple buckling model was adapted by McAdoo 
and SandwtU [1985] to the more complicated case 
for which the elastic plate possesses a brittle-ductile 
strength envelope as discussed above. They found that 
equations (1) and (2) hold quite well if the plate thick- 
ness is set equal to the thickness of the “elastic core” , 
defined as the region for which the strength of the litho- 
sphere excedes the applied stress (i.e., the thickness of 
the region between Zt and Z* in Figure 20). This ap- 
proach avoids the unreasonably high critical buckling 
stresses found when finite strength envelope is ignored. 
The critical load for this model is determined by equat- 
ing the left side of equation (1), now taken as an integral 
over the entire thickness of the crust of the supported 
stress, to the right side of equation (1), with the elas- 
tic core thickness set as the effective plate thickness. 
Since both sides of the equation depend upon the ap- 
plied stress, the critical load and corresponding buckling 
wavelength must be solved iteratively. In all of our cal- 
culations, the derived critical buckling stress is smaller 
than the maximum lithospheric strength; therefore an 
elastic core exists and buckling occurs before complete 
lithospheric failure. 

The dominant wavelength of folding determined by 
the above model depends on three parameters: the 
crustal thickness, the thermal gradient, and the strain 
rate. Figure 21 shows contour plots of the dominant 
wavelength of folding, in kilometers, as a function of 
these parameters. In each plot, one parameter is held 
constant. Figure 21a assumes a thermal gradient of 25 
K km* 1 ; Figure 216 assumes a strain rate of 10“ 16 s” 1 ; 
Figure 21c assumes a crustal thickness of 10 km. All 
results are calculated assuming a surface temperature 
of 720 K, a density for the ductile lower crust of 2800 
kg m“ 3 , a Young’s modulus of 100 GPa, and a Poisson’s 
ratio of 0.25. 

Figure 21 shows that the folding wavelength is not 
very sensitive to the crustal thickness, particularly for 
thicknesses greater than about 10 km. This is a result 
of the rapid dropoff in lithospheric strength (see Figure 
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20) due to ductile creep. The bottom of a thick crust In Figure 21 d, we have set the folding wavelength 
contributes very little strength to the crust, and there- to 8 km, a value representative of ridge spacing in the 
fore has comparatively little influence on the folding me- ridge belts of Lavinia Plan it ia. This allows us to vary 
chanics. As the crustal thickness decreases, the folding all three of the key model parameters simultaneously, 
wavelength shows more and more of a dependence, with The contours show the thermal gradient required, in 
the thinnest crust yielding the shortest wavelengths. Of kelvins per kilometer, to produce folding with a wave- 
course, the development of a second, longer wavelength length of 8 km, as a function of strain rate and crustal 


of deformation does depend strongly on crustal thick- 
ness [ Zubtr , 1987, 1990]. 

The thermal gradient and strain rate both can have a 
significant effect on the folding wavelength. Increasing 
the thermal gradient decreases the wavelength, while 
increasing the strain rate increases the wavelength. Of 
course, as Figure 21c demonstrates, the folding wave- 
length can remain essentially constant as the combina- 
tion of a relatively large strain rate and thermal gradient 
is varied to a relatively small strain rate and thermal 
gradient. This behavior is a consequence of the form 
of the creep law for the ductile deformation. Increasing 
the thermal gradient and decreasing the strain rate both 
lower the ductile yield stress of the crust, thinning the 
lithosphere and decreasing the resulting folding wave- 
length. 


thickness. Ultimately, we would like to be able to use 
this plot to constrain the thermal gradient in Lavinia 
at the time of deformation. In order to do this, how- 
ever, we need to know the strain rate involved in ridge 
belt formation. Given the geometry of the belts and 
typical estimates of the age of the Venusian surface, 
a strain rate corresponding to 10% shortening per 100 
m.y. is probably an extreme lower limit; this value is 
8 x 10~ i7 s* 1 . Taking this as the strain rate, Figure 
21 d yields a minimum thermal gradient of ~30 K kirT 1 . 
This value is somewhat greater than previous estimates 
obtained by similar arguments [Zubtr, 1987, Banerdt 
and Golombtk, 1988]. The difference is attributable 
to the ability of Magellan to resolve buckling features 
with smaller wavelengths than were inferred from low- 
resolution data. Although Figure 21 d suggests that the 





Fig. 21. Results of buckling calculations for the Venusian crust, (a) Dominant buckling wavelength as a function of 
strain rate and crustal thickness, assuming a thermal gradient of 25 K km” 1 , (b) Dominant buckling wavelength 
as a function of crustal thickness and thermal gradient, assuming a strain rate of 10“ 16 %~ l . (c) Dominant 
buckling wavelength as a function of strain rate and thermal gradient, for a crustal thickness of 10 km. (d) 
Thermal gradient as a function of strain rate and crusted thickness for a buckling wavelength (8 1cm) typical of 
the values observed in the ridge belts of Lavinia Planitia. 
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crustal thickness is not well constrained, the thickness 
of the elastic core is well constrained, with a value near 
about 200 m over the entire parameter space. 

As noted above, these results assume complete de- 
coupling between the crust and mantle. If some degree 
of coupling does exist, as we assumed in the convection 
model above, the buckling wavelengths presented here 
would increase [Banerdi and Golombck , 1988], leading 
to still higher inferred thermal gradients. Conversely, 
the buckling wavelengths would decrease if lithospheric 
strengths are reduced near the brittle-ductile transition 
by 8emibrittle deformation processes [Kirby t 1980]. Fi- 
nally, we note that there is significant uncertainty in 
the choice of creep law parameters for the Venusian 
crust and that changes in a number of parameters cho- 
sen could have a significant effect on the results. 

Conclusions 

The Magellan spacecraft has provided the first de- 
tailed look at the deformation belts of the Venusian 
lowland plains. In some ways, these observations ap- 
pear consistent with pre-Magellan expectations. The 
belts are elevated zones of intense tectonic deformation, 
and plausible arguments may be made that they are the 
consequence of crustal shortening and thickening. The 
observed scales and spacings of the belts may yield in- 
formation about the strength properties of the lower 
crust and upper mantle, and the spacings of tectonic 
features within the belts may carry information about 
the strength properties of the upper crust [Zubcr t 1987, 
1990]. However, there are some surprises as well, the 
most notable one being the fact that there are two dis- 
tinct types of belt morphology. We ultimately must un- 
derstand the reasons for these differing morphologies, 
and how they may be related to the orthogonal pat- 
tern of smaller deformational features that transects the 
plains. Perhaps the most important part of this process 
will be the examination of future Magellan data. These 
data for Lavinia will include additional SAR images at 
incidence angles that will allow stereographic determi- 
nation of topography, along with high-resolution map- 
ping of the planet’s gravity field. Moreover, images of 
other deformation belts, particularly those in Atalanta 
Planitia, have only recently become available for de- 
tailed study. A systematic comparison of the deforma- 
tion belts throughout Venus will provide a foundation 
for further modeling of belt formation and evolution. 
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